for “Animal-vectored nutrient flows over resource gradients influence the nature of local and meta-ecosystem functioning”
This notebook contains the code used to run numerical analyses of a meta-ecosystem model in which consumers move across ecosystem borders in different ways with respect to resource gradients.
In section Model Numerical Analyses, we perform numerical analyses of our model of consumer movement in a two-patch meta-ecosystem, using the set of equilibrium equations derived in Mathematica (see Mathematica notebook in this same repository; Wolfram Research Inc. (2020)). In section Changes to biomass distribution in local and meta-ecosystem, we investigate how biomass distribution changes following different types of consumer movement and how this may affect top-down influences of consumer movement on local and meta-ecosystem functions. Finally, in section Sensitivity Analyses, we perform a suite of sensitivity analyses on our results.
In this introduction, we briefly define the contrasting consumer movement types we investigate with this model (Defining movement types), discuss separation of scales in time and space (Separation of Spatial and Temporal scales), describe the state variables and parameters in the model (State variables and parameters of the model), and describe the analyses setup (Analyses setup).
Here, we focus on consumers moving over gradients of inorganic nutrients availability between two local ecosystems. Consumer movement—either gradient neutral, along-, or against-gradient—connects these two local ecosystems in a meta-ecosystem (Loreau, Mouquet, and Holt 2003). We define consumer movement as gradient neutral when it happens among ecosystems that have equal nutrient availability. Along-gradient movement takes place when consumers move following the gradient direction. Conversely, against-gradient consumer movement runs counter to the direction of the gradient. Consistent with recent reviews (Gounand et al. 2018), in the following sections we also use diffusive and non-diffusive to refer to along- and against-gradient consumer movement, respectively. In our model, we create resource gradients by varying the relative values of the inorganic nutrient input rates in the two ecosystems comprising our meta-ecosystem (parameter Ii; see Table 1 in the main text or below for a list of model parameters and state variables).
While organisms can at times cross into an ecosystem that is equally suitable for them immediately after leaving their original one, this is not often the case. Rather, organisms spend significant amounts of time searching for a new suitable ecosystem, traversing the so-called “unsuitable” matrix (Weisser and Hassell 1996). Models of organismal movement in meta-ecosystems should account for this time spent outside ecosystem borders, as it may lead to loss of biomass and nutrients from the system (Gounand et al. 2018).
We model consumer movement in our meta-ecosystem model by combining a dispersers’ pool (Weisser and Hassell 1996) with the time scales separation (TSS) mathematical technique (Otto and Day 2011). Briefly, adding a dispersers’ pool to a classic meta-ecosystem model means that consumer do not instantly travel from donor to recipient ecosystem. Instead, they move from the donor to the dispersers’ pool and then from there to the recipient (see Figure 1c in the main text for a visual representation). This allows us to separate ecosystem processes that happen locally from those that happen regionally.
Some of these ecosystem processes, like primary productivity in terrestrial ecosystems, take place over long time frames (slow processes). Others, like consumer movement, happen over short time frames (fast processes). TSS explicitly accounts for these differences by allowing fast processes to reach a quasi-equilibrium while holding slow ones constant. The quasi-equilibrium is then plugged into the equation(s) for the slow processes to investigate how these processes evolve as the fast ones stay constant (Otto and Day 2011).
The combined use of the dispersers’ pool and TSS allows us to separate local and regional dynamics in our system and to investigate the effects of consumer movement at either spatial extent at the same time. See section Model Derivation in the main text for further details on our novel approach to model consumer movement in meta-ecosystems.
The table below lists the model’s state variables and parameters, and is analogous to Table 1 in the main text. For a visual representation of the model, please see Figure 1c in the main text.
| Variable | Description | Units | Range |
|---|---|---|---|
| Ni | Nutrients stocks in patch i | g | >0 |
| Pi | Producers stocks in patch i | g | >0 |
| Ci | Consumers stocks in patch i | g | >0 |
| Q | Dispersers’ Pool | g | >0 |
| Parameter | Description | Units | Range |
|---|---|---|---|
| Ii | Inorganic nutrient input rate into patch i | g*t-1 | [0,20] |
| \(l\) | Inorganic nutrient output rate | t-1 | [0,10] |
| ui | Producer uptake rate in patch i | (g*t)-1 | [0,10] |
| ai | Consumer attack rate in patch i | (g*t)-1 | [0,10] |
| \(\epsilon\)i | Consumer assimilation efficiency in patch i | unitless | [0,1] |
| hi | Biomass loss rate from Pi | t-1 | [0,10] |
| di | Biomass loss rate from Ci | t-1 | [0,10] |
| g | Movement rate from C1 to Q | t-1 | [0,10] |
| m | Movement rate from Q to C2 | t-1 | [0,10] |
| c | Biomass loss rate from Q | t-1 | [0,10] |
We derived the model’s single feasible equilibrium in Mathematica (Wolfram Research Inc. (2020); see Model Derivation in the main text for further details).
This notebook uses the following packages:
# The following packages are used throughout this document. Before
# loading each packages, we provide aa brief description of what it
# is used for
# project dependency (packages, mostly) management
library(renv)
# data wrangling, manipulation, plotting
library(tidyverse)
# additional plotting functions
library(plotrix)
# use latin hypercube to create toy datasets
library(lhs)
# building tables in a tidyverse framework
library(gt)
# scaling functions for ggplot2
library(scales)
# color-blind palettes for ggplot2
library(khroma)
# extra tabbing functionality for rmarkdown documents
library(xaringanExtra)
# functions to time compilation of this document
library(tictoc)
# color palettes manipulation
library(RColorBrewer)
# additional color palettes inspired by MET artworks
library(MetBrewer)
# additional data wrangling
library(broom)
# print local directory structure
library(fs)
# additional ggplot2 themes
library(cowplot)
# use randomForest algorithms to perform sensitivity analyses on
# results
library(randomForest)
# extra themes for ggplot2
library(ggpubr)
# extra geoms for ggplot2
library(ggpol)
# easily create and layout multi-panel plots
library(patchwork)
# additional ggplot2 functionality
library(ggh4x)
# extra themes for ggplot2
library(ggpubr)
# extra geoms for ggplot2
library(ggpol)
# easily create and layout multi-panel plots
library(patchwork)
# Packages that need to be installed but do not require loading
# library(distill) # notebook templates library(rmarkdown) # build
# interactive documents in R
# set the theme for all ggplot2 objects
theme_set(theme_classic())
In the interest of reproducibility, we use package renv to keep track of this
notebook’s packages.
renv::snapshot()
* Lockfile written to '~/PhD/Chapter_3/Code/renv.lock'.
We compiled this notebook on a machine using the following version of R, operative system, and necessary packages:
# this needs three colons to work between package name and function call
utils:::print.sessionInfo(sessionInfo()[-8])
R version 4.2.2 (2022-10-31)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Monterey 12.6.3
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2-arm64/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices datasets utils methods
[7] base
other attached packages:
[1] ggh4x_0.2.3 patchwork_1.1.2 ggpol_0.0.7
[4] ggpubr_0.5.0 randomForest_4.7-1.1 cowplot_1.1.1
[7] fs_1.6.0 broom_1.0.3 MetBrewer_0.2.0
[10] RColorBrewer_1.1-3 tictoc_1.1 xaringanExtra_0.7.0
[13] khroma_1.9.0 scales_1.2.1 gt_0.8.0
[16] lhs_1.1.6 plotrix_3.8-2 forcats_1.0.0
[19] stringr_1.5.0 dplyr_1.1.0 purrr_1.0.1
[22] readr_2.1.3 tidyr_1.3.0 tibble_3.1.8
[25] ggplot2_3.4.0 tidyverse_1.3.2 renv_0.16.0
Finally, this notebook assumes the following folder structure:
dir_tree(path = "../", recurse = 0, all = FALSE, type = "directory")
../
├── Code
├── Data
├── LitRev
├── Manuscripts
├── Modeling_Ch3
├── Results
└── Submission
This structure allows us to take advantage of relative paths when loading data
or saving figures. Note that ../ stands for the root folder. If you wish to
recompile this notebook, please organize your folder structure
accordingly.
Here, we perform the analyses of the meta-ecosystem model described in section Numerical Analyses of the main text. If you are interested in the code used to produce the main text’s figures and tables, this can be found in sections Visual Deliverables—Ecosystem functions and Summary Tables—Ecosystem functions below.
We conduct three rounds of analyses, one each for the three types of consumer movement we are interested in: gradient-neutral, along-, and against-gradient. Each round of analysis comprises iterating the model over n = 10 000 parameter sets, running stability analyses, and graphing the results. The analyses will follow these steps:
This notebooks defaults to loading the two pre-generated datasets provided
alongside it, in folder ../Data. DATA1 includes parameters scaled [0–10],
whereas DATA2 includes parameters scaled [0–1]. These were generated with a
Latin hypercube sampling design, using function lhs(), and then used to
produce the results and figures presented in the paper. To generate new randomly
drawn parameter values, set eval=FALSE in the data-load code chunk below
and eval=TRUE in the following one, labelled disp-poolsetParams.
In the code chunk below, we show how to generate two sets of randomly-drawn
parameter values to use in all scenarios. We use function lhs() to do so. One
set includes parameters values scaled 0–10 (DATA1), while the other one
includes parameter values scaled 0–1 (DATA2).
# Below we generate random values for those parameters that will vary
# during the simulations: attack rates, efficiency, death rates, and
# movement rates
DATA1 <- NULL
# Run 2000 iterations of 100 numbers (in 100 bins between 0 to 1)
for (i in seq(1, 100, 1)) {
# Use lhs function in R. 100 bins, 11 iterations - one for each
# parameter that is scaled between 0 and 10.
loop_1 <- randomLHS(100, 11)
# Use rbind to stack each iteration on top of the other to create
# a one column list of 10000 samples with equal distribution in
# 100 bins. Multiplied by 10 to get range of data from 0 to 10
DATA1 <- rbind(DATA1, data.frame(i = i, Result = loop_1 * 10))
# print the i value so you can track progress print(i) Close the
# loop
}
# This is for 2 parameters that is scaled 0 to 1.
DATA2 <- NULL
# Run 2000 iterations of 100 numbers (in 100 bins between 0 to 1)
for (i in seq(1, 100, 1)) {
# Use lhs function in R. 100 bins, 2 iterations.
loop_2 <- randomLHS(100, 2)
# Use rbind to stack each iteration on top of the other to create
# a one column list of 200000 samples with equal distribution in
# 100 bins. Multiplied by 10 to get range of data from 0 to 10
DATA2 <- rbind(DATA2, data.frame(i = i, Result = loop_2))
# print the i value so you can track progress print(i) Close the
# loop
}
# Save these two datasets for later re-use write.csv(DATA1, file =
# '../Data/DATA1.csv', row.names = FALSE) write.csv(DATA2, file =
# '../Data/DATA2.csv', row.names = FALSE)
We begin our numerical analyses of our model with the control scenario, in which consumers move between two ecosystems that have equal environmental fertility. Consumer movement in this scenario is gradient-neutral, as the environment is homogeneous in terms of resource availability. The results of this first set of model iterations produce the values used at the denominator of the response ratio (LRR) used to capture changes in local and meta-ecosystem function values (see section [Visual Deliverables] below). Below, we describe each step of the analyses in detail: later analyses will not describe these steps again.
First, we allocate an empty dataframe PRED to contain the results of each
iteration of the model. This dataframe will contain the parameter values and the
resulting local and meta-ecosystem function values.
# Allocate an empty data frame to save the data in, then allocate
# some rows and columns. We will replace the current numbers with
# results below Data we are calculating are (bio)mass for soil
# nutrients (N1, N2), producers (P1, P2), consumers (C1, C2),
# disperses (Q)
PRED <- NULL
PRED <- rbind(PRED, data.frame(TIME = seq(0, 100, 1), I1 = seq(0, 100,
1), I2 = seq(0, 100, 1), l = seq(0, 100, 1), u1 = seq(0, 100, 1), u2 = seq(0,
100, 1), a1 = seq(0, 100, 1), a2 = seq(0, 100, 1), h1 = seq(0, 100,
1), h2 = seq(0, 100, 1), d1 = seq(0, 100, 1), d2 = seq(0, 100, 1),
g = seq(0, 100, 1), m = seq(0, 100, 1), c = seq(0, 100, 1), e1 = seq(0,
100, 1), e2 = seq(0, 100, 1), N1 = seq(0, 100, 1), P1 = seq(0,
100, 1), C1 = seq(0, 100, 1), N2 = seq(0, 100, 1), P2 = seq(0,
100, 1), C2 = seq(0, 100, 1), Q = seq(0, 100, 1), FLUX_P1 = seq(0,
100, 1), FLUX_C1 = seq(0, 100, 1), FLUX_P2 = seq(0, 100, 1), FLUX_C2 = seq(0,
100, 1), FLUX_Eco_1 = seq(0, 100, 1), FLUX_Eco_2 = seq(0, 100,
1), FLUX_Q = seq(0, 100, 1), PROD_P1 = seq(0, 100, 1), PROD_C1 = seq(0,
100, 1), PROD_P2 = seq(0, 100, 1), PROD_C2 = seq(0, 100, 1)))
Now we run the model. We first set the values for both nutrient input rates in
ecosystem 1 (I1 in the code chunks) and ecosystem 2 (I2) at 10. Recall that
this is the control scenario with equal environmental fertility. We also set the
value of the inorganic nutrient leaching rate (l) at 0.1.
The code chunk below runs the model using a for loop that, over 10 000 time
steps, iterates the model’s single feasible equilibrium to calculate:
N1, N2), Producers (P1, P2), and Consumers (C1,
C2)Q) stockFLUX_X, where X is the compartment)PROD_X, where X is the
trophic compartment)in both local ecosystems. At each time step, we draw a random value for each
parameter from either the DATA1 or DATA2 dataset—depending on whether the
parameter is scaled 0–10 or 0–1, respectively. Then, we plug these values into
the model’s equilibria and ecosystem function formulas derived in Mathematica
(Wolfram Research Inc. 2020), and finally store both the randomly drawn parameter values and
ecosystem function values in PRED.
Note that feasibility conditions for this model cannot be resolved analytically. Hence, we do not check that our random parameters values satisfy them. Rather, we run stability analyses of each parameter set below.
# Here we set values for parameters that we will not vary during the simulations
# Nutrients input rate in each patch
I1 = 10
I2 = 10
# Patch leaching rate/soil loss rate
l = 0.1
# We want to start simulations at 0
i1 = 0
for (i in seq(0,9999,1)){
# go to the next row
i1 = i1 + 1
# Let's sample the parameter values from DATA. We will sample these i times as
# defined by the for loop above
# producers uptake rates
u1 = DATA1[i1,2] # in ecosystem 1
u2 = DATA1[i1,3] # in ecosystem 2
# consumers attack rates
a1 = DATA1[i1,4] # in ecosystem 1
a2 = DATA1[i1,5] # in ecosystem 2
# death rates (= recycling rate)
h1 = DATA1[i1,6] # death (=recycling) rate for producers in ecosystems 1
h2 = DATA1[i1,7] # death (=recycling) rate for producers in ecosystems 2
d1 = DATA1[i1,8] # death (=recycling) rate for consumers in ecosystem 1
d2 = DATA1[i1,9] # death (=recycling) rate for consumers in ecosystem 2
c = DATA1[i1,10] # death rate in the disperser's pool Q
# movement rate to the Dispersers' pool Q
g = DATA1[i1,11]
# movement rate from the Dispersers' pool Q
m = DATA1[i1,12]
# assimilation efficiency of consumers
e1 = DATA2[i1,2] # in ecosystem 1
e2 = DATA2[i1,3] # in ecosystem 2
PRED[i1,] <- c(i1, I1, I2, l, u1, u2, a1, a2, h1, h2, d1, d2, g, m, c, e1, e2,
# Nutrient stock in ecosystem 1
((d1 - d1*e1 + g)*h1 + a1*e1*I1)/(a1*e1*l + (d1 - d1*e1 + g)*u1),
# Primary Producers biomass in ecosystem 1
(d1 + g)/(a1*e1),
# Consumers biomass in ecosystem 1
(e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
# Nutrient stock in ecosystem 2
(-a1*e1*(d2*(-1 + e2)*h2 - a2*e2*I2)*l*(c + m) + d2*(-1 + e2)*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*(e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1))))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2)),
# Primary Producers biomass in ecosystem 2
(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2),
# Consumer biomass in ecosystem 2
((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2))/(a2*e2*l - d2*(-1 + e2)*u2),
# Consumer biomass in the Dispersers' Pool (Q)
(e1*g*(-h1*l + I1*u1))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
# Nutrient Flux for Primary Producers in Ecosystem 1
((d1 + g)*h1)/(a1*e1),
# Nutrient Flux for Consumers in Ecosystem 1
(d1*e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
# Nutrient Flux for Primary Producers in Ecosystem 2
(h2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2),
# Nutrient Flux for Consumers in Ecosystem 2
(d2*((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2)))/(a2*e2*l - d2*(-1 + e2)*u2),
# Nutrient Flux for Ecosystem 1
((d1 + g)*h1)/(a1*e1) + (d1*e1*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
# Nutrient Flux for Ecosystem 2
(h2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2) + (d2*((e1*g*m*(-h1*l + I1*u1)*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)) + e2*(-h2*l + I2*u2)))/(a2*e2*l - d2*(-1 + e2)*u2),
# Nutrient Flux in the Dispersers' Pool (Q)
(c*e1*g*(-h1*l + I1*u1))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
# Primary Productivity in Ecosystem 1
((d1 + g)*((d1 - d1*e1 + g)*h1 + a1*e1*I1)*u1)/(a1*e1*(a1*e1*l + (d1 - d1*e1 + g)*u1)),
# Consumer (=Secondary) Productivity in Ecosystem 1
(e1*(d1 + g)*(-h1*l + I1*u1))/(a1*e1*l + (d1 - d1*e1 + g)*u1),
# Primary Productivity in Ecosystem 2
((-a1*e1*(d2*(-1 + e2)*h2 - a2*e2*I2)*l*(c + m) + d2*(-1 + e2)*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*(e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1))))*u2*(l*(-d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(-h1*l + I1*u1)) + d2*(d1*e1*I2*(c + m)*u1 - (d1 + g)*I2*(c + m)*u1 + e1*g*m*(h1*l - I1*u1))*u2 + a1*d2*e1*l*(c + m)*(h2*l - I2*u2)))/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2)*(a2*e2*h2*l*(c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1) - a2*(a1*e1*e2*I2*l*(c + m) + e2*(d1 + g)*I2*(c + m)*u1 - e1*(d1*e2*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2)),
# Consumer (=Secondary) Productivity in Ecosystem 2
(e2*l*(-a1*d2*e1*h2*l*(c + m) + d2*(d1*(-1 + e1) - g)*h2*(c + m)*u1 + a2*e1*g*m*(h1*l - I1*u1)) + d2*e2*(a1*e1*I2*l*(c + m) + (d1 + g)*I2*(c + m)*u1 - e1*(d1*I2*(c + m)*u1 + g*m*(h1*l - I1*u1)))*u2)/((c + m)*(a1*e1*l + (d1 - d1*e1 + g)*u1)*(a2*e2*l - d2*(-1 + e2)*u2))
)
}
To investigate any influence of consumer movement on nutrient cycling
at regional spatial scales, we calculate local and meta-ecosystem total biomass
(Ntot,Ptot, Ctot), meta-ecosystem nutrient flux (MetaEcoFlux) and
trophic compartment nutrient flux (FLUX_Ptot, FLUX_Ctot) and productivity
(PROD_Ptot, PROD_Ctot). Furthermore, we manually calculate each ecosystem’s
own total recycling flux as the sum of the flux of its 2 biotic compartments, as
a check on the Mathematica-derived formulas used above.
# calculate meta-ecosystem biomass stocks for each trophic
# compartment, a.k.a., meta-ecosystem biomass
PRED$Ntot <- PRED$N1 + PRED$N2
PRED$Ptot <- PRED$P1 + PRED$P2
PRED$Ctot <- PRED$C1 + PRED$C2
# now calculate recycling flux for the local ecosystem manually to
# double check the formulae used above are correct
PRED$FLUX_Eco1_check <- PRED$FLUX_P1 + PRED$FLUX_C1
PRED$FLUX_Eco2_check <- PRED$FLUX_P2 + PRED$FLUX_C2
# and the meta-ecosystem recycling flux
PRED$FLUX_Ptot <- PRED$FLUX_P1 + PRED$FLUX_P2
PRED$FLUX_Ctot <- PRED$FLUX_C1 + PRED$FLUX_C2
PRED$MetaEcoFlux <- PRED$FLUX_Eco_1 + PRED$FLUX_Eco_2 - PRED$FLUX_Q
# finally, calculate the production for each compartment at
# meta-ecosystem level
PRED <- dplyr::mutate(PRED, PROD_Ptot = PROD_P1 + PROD_P2, .after = "PROD_C2")
PRED <- dplyr::mutate(PRED, PROD_Ctot = PROD_C1 + PROD_C2, .after = "PROD_Ptot")
Parameter sets that produce values of the model’s state variable that are \(\leq 0\) lack biological meaning. Accordingly, in the code chunk below, we flag these parameter sets for later exclusion from the analysis.
# Equilibrium stocks less than or equal to 0 make no biological
# sense, so we will flag them in the PRED dataset.
PRED <- PRED %>%
dplyr::mutate(., biosense = ifelse(N1 > 0 & P1 > 0 & C1 > 0 & N2 >
0 & P2 > 0 & C2 > 0 & Q > 0, "yes", "no"), .after = MetaEcoFlux) %>%
dplyr::mutate_at(vars(biosense), factor)
We perform the stability analyses on the model by fitting the same sets of parameter values and equilibrium state variable values used above to the partial derivatives extracted from the model’s Jacobian matrix evaluated in Mathematica (Wolfram Research Inc. 2020). We evaluate stability of each parameter set by checking whether the real part of the leading eigenvalue of the Jacobian is positive or negative (Otto and Day 2011). If the real part of the Jacobian’s leading eigenvalue is negative, the equilibrium is stable (Otto and Day 2011).
# create an empty dataframe
MathStab <- NULL
for (i in 1:nrow(PRED)) {
# fetch parameter values from the equilibrium simulations df
I1 = PRED$I1[i]
I2 = PRED$I2[i]
l = PRED$l[i]
u1 = PRED$u1[i]
u2 = PRED$u2[i]
a1 = PRED$a1[i]
a2 = PRED$a2[i]
h1 = PRED$h1[i]
h2 = PRED$h2[i]
d1 = PRED$d1[i]
d2 = PRED$d2[i]
g = PRED$g[i]
m = PRED$m[i]
c = PRED$c[i]
e1 = PRED$e1[i]
e2 = PRED$e2[i]
# Put in the solutions to the equilibrium values here:
dN1 = PRED$N1[i]
dP1 = PRED$P1[i]
dC1 = PRED$C1[i]
dN2 = PRED$N2[i]
dP2 = PRED$P2[i]
dC2 = PRED$C2[i]
# create the Jacobian from Mathematica NOTE: in transposing from
# Mathematica, the Jacobian is here assembled by row
Jacob <- rbind(c(-l - dP1 * u1, h1 - dN1 * u1, d1, 0, 0, 0), c(dP1 *
u1, -a1 * dC1 - h1 + dN1 * u1, -a1 * dP1, 0, 0, 0), c(0, a1 * dC1 *
e1, -d1 - g + a1 * e1 * dP1, 0, 0, 0), c(0, 0, 0, -l - dP2 * u2,
h2 - dN2 * u2, d2), c(0, 0, 0, dP2 * u2, -a2 * dC2 - h2 + dN2 *
u2, -a2 * dP2), c(0, 0, (g * m)/(c + m), 0, a2 * dC2 * e2, -d2 +
a2 * e2 * dP2))
# Check to see if the real part of the leading eigenvalue is less
# than 0 or not - if it is than the system is stable.
if (max(Re(base::eigen(Jacob)$values)) < 0) {
stable <- "stable"
} else {
stable <- "unstable"
}
MathStab <- rbind(MathStab, data.frame(TIME = i, I1, I2, l, u1, u2,
a1, a2, h1, h2, d1, d2, g, m, c, e1, e2, dN1, dP1, dC1, dN2, dP2,
dC2, EiV1 = Re(eigen(Jacob)$values[1]), EiV2 = Re(eigen(Jacob)$values[2]),
EiV3 = Re(eigen(Jacob)$values[3]), EiV4 = Re(eigen(Jacob)$values[4]),
EiV5 = Re(eigen(Jacob)$values[5]), EiV6 = Re(eigen(Jacob)$values[6]),
maxEv = max(Re(base::eigen(Jacob)$values)), stable = stable, biosense = PRED$biosense[i]))
}
MathStab$stable <- as.factor(MathStab$stable)
# separate unstable equilibria to work with later
MathStabUS <- subset(MathStab, MathStab$stable == "unstable")
According to the stability analyses run above, about 1.55% of the parameter sets produced unstable results (i.e., 155 out of 10000 iterations).
Now, let’s check if the parameter sets that produce unstable equilibria are also those that produce equilibria lacking biological meaning—i.e., those where state variables values at equilibrium are \(\leq 0\).
biononsense <- subset(PRED, PRED$biosense == "no")
identical(as.numeric(MathStabUS[, "TIME"]), biononsense[, "TIME"])
[1] TRUE
It appears that they are. Hence, let’s add the stable/unstable information to
PRED.
Now, let’s exclude these unstable parameter sets from the analyses below and from our results.
We are also going to sample the dataset
PRED, that contains the results of our model’s iterations, for 1000 random
iterations results. We will store these in PRED_1k and use them in Figures
1, 2, 3,
4, 5 below
to produce some preliminary graphs of our results.
PREDpos <- filter(PRED, PRED$biosense == "yes" & PRED$stable == "stable")
# sample PREDpos to only use 1000 random simulation results
PRED_sample1000 <- PREDpos[sample(nrow(PREDpos), size = 1000), ]
PRED_1k <- droplevels(PRED_sample1000)
Before proceeding further, let’s save the PRED object to disk, for ease of
use in the future.
write.csv(PRED, file = "../Results/PRED.csv", row.names = FALSE)
We now investigate how the presence of unidirectional, gradient-neutral consumer
movement influences the meta-ecosystem’s behavior using box-plots produced using
PRED_1k. Note that in the figures below, “Donor” always indicates ecosystem 1
and “Recipient” always indicates ecosystem 2.
# pivot the dataset to longer format for boxplot plotting
PRED_biomass_long <- select(PRED_1k, N1:Q, Ntot:Ctot) %>%
pivot_longer(., names_to = "Compartment", values_to = "Stock", cols = c(N1,
P1, C1, N2, P2, C2, Q, Ntot, Ptot, Ctot)) %>%
dplyr::mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
"C1" | Compartment == "P1", "Donor", if_else(Compartment == "N2" |
Compartment == "P2" | Compartment == "C2", "Recipient", if_else(Compartment ==
"Q", "Dispersers Pool", "Meta-ecosystem"))))
PRED_biomass_long$Compartment <- as_factor(PRED_biomass_long$Compartment)
PRED_biomass_long$Ecosystem <- as_factor(PRED_biomass_long$Ecosystem)
comparts_medians <- PRED_biomass_long %>%
dplyr::group_by(Compartment, Ecosystem) %>%
dplyr::summarise(., median = median(Stock), .groups = "keep")
ggplot(data = PRED_biomass_long, aes(x = Compartment, y = Stock, col = Ecosystem)) +
geom_boxjitter(alpha = 0.25, outlier.intersect = TRUE, outlier.shape = 25,
jitter.shape = 21, outlier.stroke = 1, jitter.stroke = 1) + stat_summary(geom = "label",
fun = quantile, fun.args = list(probs = 0.5), aes(label = round(after_stat(y),
2)), position = position_dodge(width = 0.2), show.legend = F, size = 3,
label.padding = unit(0.15, "lines"), label.size = 0.15, fontface = "bold",
color = "black") + facet_grid(. ~ Ecosystem, scales = "free") + scale_color_manual(values = met.brewer("Egypt",
4)) + theme_pubr() + coord_cartesian(y = c(0, 100)) + theme(legend.position = "none")
Figure 1: Biomass stock values at equilibrium for local and meta-ecosystem scale, drawn using 1k randomly-drawn data points from the model simulations. Consumer movement from ecosystem 1 to ecosystem 2 via the dispersers’ pool results in foraging pressure of consumers being transferred from producers in ecosystem 1 to those in ecosystem 2. Primary producers in ecosystem 1 are thus released from foraging pressure while those in ecosystem 2 are depressed. In turn, this leads to differences in nutrient stocks between the two ecosystems at both the local and meta-ecosystem extents. Point-down triangles identify outliers; labels on top of the compartments’ boxplots report median values for the full set of simulations (i.e., 10k iterations). Values on the y-axis are limited between [0, 100] to zoom in on the variation in the results.
# ggsave(filename = '../Results/ModelEq_base.pdf', device = 'pdf',
# dpi = 300, height = 5)
The graph above shows a spatial trophic cascade (Knight et al. 2005). In ecosystem 2,
the influx of consumers from ecosystem 1 increases the consumer biomass, in turn
depressing the producers’ abundance and releasing the nutrient stock from the
producers’ trophic pressure. Conversely, in ecosystem 1, consumers’ numbers fall
steadily—due to the combined action of mortality and emigration towards
ecosystem 2—thus releasing the producers from the trophic pressure exerted by
consumers. In turn, this lets producers grow unchecked, depressing the nutrient
stocks in ecosystem 1. We can see this as well by looking at the change of log10
biomass values for each trophic compartment with respect to the two movement
parameters: g, the rate of consumer movement from ecosystem 1 to the
dispersers’ pool Q (Figure 2), and m, the rate of consumer
movement from the dispersers’ pool Q to ecosystem 2 (Figure 3).
eco_levels <- c(N1 = "Nutrient", P1 = "Producers", C1 = "Consumers", N2 = "Nutrient",
P2 = "Producers", C2 = "Consumers", Q = "Dispersers", Ntot = "Nutrient",
Ptot = "Producers", Ctot = "Consumers")
# pivot the dataset to longer format for boxplot plotting
PRED_long <- pivot_longer(PRED_1k, names_to = "Compartment", values_to = "Stock",
cols = c(N1, P1, C1, N2, P2, C2, Q, Ntot, Ptot, Ctot)) %>%
dplyr::mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
"C1" | Compartment == "P1", "Donor", if_else(Compartment == "N2" |
Compartment == "P2" | Compartment == "C2", "Recipient", if_else(Compartment ==
"Q", "Dispersers", "Meta-ecosystem"))))
PRED_long$Compartment <- as_factor(PRED_long$Compartment)
PRED_long$Ecosystem <- as_factor(PRED_long$Ecosystem)
ggplot(PRED_long, aes(g, log10(Stock), col = Ecosystem)) + geom_point(alpha = 0.15,
shape = 21, size = 2, stroke = 1) + geom_smooth(method = "lm", col = "white",
se = T, lwd = 0.4) + facet_grid(. ~ Compartment) + scale_color_manual(values = met.brewer("Egypt",
4)) + theme_pubr() + guides(colour = guide_legend(override.aes = list(alpha = 1))) +
theme(axis.text.x = element_text(angle = 45, vjust = 0.5, size = 8),
text = element_text(size = 11))
Figure 2: Change in biomass in the meta-ecosystems’ trophic compartments with change in the movement rate from the donor ecosystem 1 to the dispersers’ pool (g). Each facet presents results for a trophic compartment, with color identifying local and meta-ecosystem. White lines are regression lines with 95% CI (grey shaded areas). Note the log10 scale on the y-axis. Graph produced using 1000 randomly-drawn data points from the model’s analyses.
ggplot(PRED_long, aes(m, log10(Stock), col = Ecosystem)) + geom_point(alpha = 0.15,
shape = 21, size = 2, stroke = 1) + geom_smooth(method = "lm", col = "white",
se = T, lwd = 0.4) + facet_grid(. ~ Compartment) + scale_color_manual(values = met.brewer("Egypt",
4)) + theme_pubr() + guides(colour = guide_legend(override.aes = list(alpha = 1))) +
theme(axis.text.x = element_text(angle = 45, vjust = 0.5, size = 8),
text = element_text(size = 11))
Figure 3: Change in biomass in the meta-ecostystem’s trophic compartments with change in movement rate from the disperser pool to the recipient ecosystem 2 (m). Graph produced using 1000 randomly-drawn data points from the model’s analyses. Note the log10 scale on the y-axis. All other specifications as in Figure 2.
Finally, we can see a similar effect when looking at the change in biomass in
the two local ecosystems with change in the mortality rate of consumers while in
the dispersers’ pool (parameter c in the model; Figure 4).
ggplot(PRED_long, aes(c, log10(Stock), col = Ecosystem)) + geom_point(alpha = 0.15,
shape = 21, size = 2, stroke = 1) + geom_smooth(method = "lm", col = "white",
se = T, lwd = 0.4) + facet_grid(. ~ Compartment) + scale_color_manual(values = met.brewer("Egypt",
4)) + theme_pubr() + guides(colour = guide_legend(override.aes = list(alpha = 1))) +
theme(axis.text.x = element_text(angle = 45, vjust = 0.5, size = 8),
text = element_text(size = 11))
Figure 4: Change in biomass in the meta-ecostystem’s trophic compartments with change in the consumer death rate while in the dispersers’ pool (c). This could represent predation risk, but also environmental hazards and stochastic events such as crossing a dangerous linear feature of the landscape—e.g., a busy highway. Graph produced using 1000 randomly-drawn data points from the model’s analyses. Note the log10 scale on the y-axis. All other specifications as in Figure 2.
In the graphs below (Figure 5), we look at nutrient flux values in each ecosystem and in the meta-ecosystem for each iteration of the model.
First, we will extract the flux columns (see above) from the PRED_1k dataframe
into a new one which we will call FluxPRED. We will then pivot this dataframe
from wide to long format, using function pivot_longer from package tidyr.
FluxPRED <- PRED_1k %>%
dplyr::select(TIME:c, FLUX_P1:FLUX_Eco_2, MetaEcoFlux) %>%
tidyr::pivot_longer(., names_to = "Compartment", values_to = "Flux",
cols = c(FLUX_P1:MetaEcoFlux)) %>%
dplyr::mutate(., Scale = if_else(Compartment == "FLUX_P1" | Compartment ==
"FLUX_C1" | Compartment == "FLUX_Eco_1", "Donor", if_else(Compartment ==
"FLUX_P2" | Compartment == "FLUX_C2" | Compartment == "FLUX_Eco_2",
"Recipient", "Meta-ecosystem"))) %>%
dplyr::mutate(., Scale = fct_relevel(Scale, "Donor", "Recipient", "Meta-ecosystem"))
Now, we produce the box-and-whiskers plot of recycling flux in local and meta-ecosystem.
FluxPRED$Scale <- as.factor(FluxPRED$Scale)
FluxPRED$Compartment <- as.factor(FluxPRED$Compartment)
comparts_medians <- FluxPRED %>%
dplyr::group_by(Compartment, Scale) %>%
dplyr::summarise(., median = median(Flux), .groups = "keep")
FluxPRED %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, c("FLUX_P1",
"FLUX_C1", "FLUX_Eco_1", "FLUX_P2", "FLUX_C2", "FLUX_Eco_2", "MetaEcoFlux"))) %>%
dplyr::mutate(., Scale = fct_relevel(Scale, c("Donor", "Recipient",
"Meta-ecosystem"))) %>%
ggplot(., aes(x = Compartment, y = Flux)) + geom_boxjitter(outlier.intersect = TRUE,
alpha = 0.25, outlier.shape = 25, jitter.shape = 21) + stat_summary(geom = "label",
fun = quantile, fun.args = list(probs = 0.5), aes(label = round(after_stat(y),
2)), position = position_dodge(width = 0.2), show.legend = F, size = 3,
label.padding = unit(0.15, "lines"), label.size = 0.15, fontface = "bold",
color = "black") + facet_grid(. ~ Scale, scales = "free") + scale_x_discrete(labels = c(FLUX_C1 = "C1",
FLUX_P1 = "P1", FLUX_Ptot = "Producers", FLUX_C2 = "C2", FLUX_P2 = "P2",
FLUX_Ctot = "Consumers", FLUX_Eco_1 = "Donor", FLUX_Eco_2 = "Recipient",
MetaEcoFlux = "Meta-ecosystem")) + theme_pubr() + coord_cartesian(y = c(0,
500)) + theme(legend.position = "none", axis.title.x = element_blank())
Figure 5: Nutrient flux at local and meta-ecosystem scales when consumers move from ecosystem 1 to 2. The two ecosystems have the same fertility—i.e., basal inorganic input values are the same. Point-down triangles identify outliers; labels on top of the compartments’ boxplots report median values for the full set of simulations. The y-axis is constrained between [0, 500] to show the variation in the results. Graph produced using 1000 randomly-drawn data points from the model’s analyses.
Here, we investigate the influence of consumers moving from ecosystem 1 to
ecosystem 2 when there is a gradient of nutrient availability between the two
ecosystems. We will investigate two scenarios: higher environmental fertility in
the donor ecosystem 1 (i.e., I1 >> I2) and higher environmental fertility in
the recipient ecosystem 2 (i.e., I1 << I2). Note that the environmental
leaching rate of each ecosystem—i.e., the rate at which nutrients
leave local ecosystems independently of trophic processes—is invariant and equal
between the two ecosystems (l = 0.1).
The analysis will follow the same steps as before: draw parameter values
from our DATA1 and DATA2 data objects, simulate model
behavior, run stability analyses, remove equilibria that are unstable and/or
lack biological sense, and graphically inspect and interpret the results.
We begin by simulating a meta-ecosystem where the donor ecosystem 1 has higher environmental fertility than the recipient ecosystem 2. Consumers move along-gradient, in a diffusive way that resembles how organismal movement has been represented in previous studies (e.g., Gravel et al. 2010).
# Here we set values for inorganic Nutrients input rate in each patch
# In this case, I1 > I2
I1 = 18
I2 = 2
# Leaching does not change as is still equal across ecosystems
l = 0.1
# Allocate an empty data frame to save the data in
PREDI2 <- NULL
# Allocate some rows and columns. We will replace the current numbers
# with results below Data we are calculating are (bio)mass for soil
# nutrients (N1, N2), producers (P1, P2), consumers (C1,C2),
# disperses (Q)
PREDI2 <- rbind(PREDI2, data.frame(TIME = seq(0, 100, 1), I1 = I1, I2 = I2,
l = l, u1 = seq(0, 100, 1), u2 = seq(0, 100, 1), a1 = seq(0, 100, 1),
a2 = seq(0, 100, 1), h1 = seq(0, 100, 1), h2 = seq(0, 100, 1), d1 = seq(0,
100, 1), d2 = seq(0, 100, 1), g = seq(0, 100, 1), m = seq(0, 100,
1), c = seq(0, 100, 1), e1 = seq(0, 100, 1), e2 = seq(0, 100, 1),
N1 = seq(0, 100, 1), P1 = seq(0, 100, 1), C1 = seq(0, 100, 1), N2 = seq(0,
100, 1), P2 = seq(0, 100, 1), C2 = seq(0, 100, 1), Q = seq(0, 100,
1), FLUX_P1 = seq(0, 100, 1), FLUX_C1 = seq(0, 100, 1), FLUX_P2 = seq(0,
100, 1), FLUX_C2 = seq(0, 100, 1), FLUX_Eco_1 = seq(0, 100, 1),
FLUX_Eco_2 = seq(0, 100, 1), FLUX_Q = seq(0, 100, 1), PROD_P1 = seq(0,
100, 1), PROD_C1 = seq(0, 100, 1), PROD_P2 = seq(0, 100, 1), PROD_C2 = seq(0,
100, 1)))
# We want to start simulations at 0
i1 = 0
for (i in seq(0, 9999, 1)) {
# go to the next row
i1 = i1 + 1
# Let's sample the parameter values from DATA. We will sample
# these x times as defined by the i for loop above producers
# uptake rates
u1 = DATA1[i1, 2] # in ecosystem 1
u2 = DATA1[i1, 3] # in ecosystem 2
# consumers attack rates
a1 = DATA1[i1, 4]
a2 = DATA1[i1, 5]
# death rates (= recycling rate)
h1 = DATA1[i1, 6] # death (=recycling) rate for producers in ecosystems 1
h2 = DATA1[i1, 7] # death (=recycling) rate for producers in ecosystem 2
d1 = DATA1[i1, 8] # death (=recycling) rate for consumers in ecosystem 1
d2 = DATA1[i1, 9] # death (=recycling) rate for consumers in ecosystem 2
c = DATA1[i1, 10] # death rate in the disperser's pool Q
# movement rate to the Disperser's pool Q
g = DATA1[i1, 11]
# movement rate from the Disperser's pool Q
m = DATA1[i1, 12]
# efficiency of consumers
e1 = DATA2[i1, 2] # in ecosystem 1
e2 = DATA2[i1, 3] # in ecosystem 2
PREDI2[i1, ] <- c(i1, I1, I2, l, u1, u2, a1, a2, h1, h2, d1, d2, g,
m, c, e1, e2, ((d1 - d1 * e1 + g) * h1 + a1 * e1 * I1)/(a1 * e1 *
l + (d1 - d1 * e1 + g) * u1), (d1 + g)/(a1 * e1), (e1 * (-h1 *
l + I1 * u1))/(a1 * e1 * l + (d1 - d1 * e1 + g) * u1), (-a1 *
e1 * (d2 * (-1 + e2) * h2 - a2 * e2 * I2) * l * (c + m) + d2 *
(-1 + e2) * (d1 * (-1 + e1) - g) * h2 * (c + m) * u1 + a2 *
(e2 * (d1 + g) * I2 * (c + m) * u1 - e1 * (d1 * e2 * I2 * (c +
m) * u1 + g * m * (h1 * l - I1 * u1))))/((c + m) * (a1 *
e1 * l + (d1 - d1 * e1 + g) * u1) * (a2 * e2 * l - d2 * (-1 +
e2) * u2)), (l * (-d2 * (d1 * (-1 + e1) - g) * h2 * (c + m) *
u1 + a2 * e1 * g * m * (-h1 * l + I1 * u1)) + d2 * (d1 * e1 *
I2 * (c + m) * u1 - (d1 + g) * I2 * (c + m) * u1 + e1 * g *
m * (h1 * l - I1 * u1)) * u2 + a1 * d2 * e1 * l * (c + m) *
(h2 * l - I2 * u2))/(a2 * e2 * h2 * l * (c + m) * (a1 * e1 *
l + (d1 - d1 * e1 + g) * u1) - a2 * (a1 * e1 * e2 * I2 * l *
(c + m) + e2 * (d1 + g) * I2 * (c + m) * u1 - e1 * (d1 * e2 *
I2 * (c + m) * u1 + g * m * (h1 * l - I1 * u1))) * u2), ((e1 *
g * m * (-h1 * l + I1 * u1) * u2)/((c + m) * (a1 * e1 * l +
(d1 - d1 * e1 + g) * u1)) + e2 * (-h2 * l + I2 * u2))/(a2 *
e2 * l - d2 * (-1 + e2) * u2), (e1 * g * (-h1 * l + I1 * u1))/((c +
m) * (a1 * e1 * l + (d1 - d1 * e1 + g) * u1)), ((d1 + g) *
h1)/(a1 * e1), (d1 * e1 * (-h1 * l + I1 * u1))/(a1 * e1 * l +
(d1 - d1 * e1 + g) * u1), (h2 * (l * (-d2 * (d1 * (-1 + e1) -
g) * h2 * (c + m) * u1 + a2 * e1 * g * m * (-h1 * l + I1 *
u1)) + d2 * (d1 * e1 * I2 * (c + m) * u1 - (d1 + g) * I2 *
(c + m) * u1 + e1 * g * m * (h1 * l - I1 * u1)) * u2 + a1 *
d2 * e1 * l * (c + m) * (h2 * l - I2 * u2)))/(a2 * e2 * h2 *
l * (c + m) * (a1 * e1 * l + (d1 - d1 * e1 + g) * u1) - a2 *
(a1 * e1 * e2 * I2 * l * (c + m) + e2 * (d1 + g) * I2 * (c +
m) * u1 - e1 * (d1 * e2 * I2 * (c + m) * u1 + g * m * (h1 *
l - I1 * u1))) * u2), (d2 * ((e1 * g * m * (-h1 * l + I1 *
u1) * u2)/((c + m) * (a1 * e1 * l + (d1 - d1 * e1 + g) * u1)) +
e2 * (-h2 * l + I2 * u2)))/(a2 * e2 * l - d2 * (-1 + e2) *
u2), ((d1 + g) * h1)/(a1 * e1) + (d1 * e1 * (-h1 * l + I1 *
u1))/(a1 * e1 * l + (d1 - d1 * e1 + g) * u1), (h2 * (l * (-d2 *
(d1 * (-1 + e1) - g) * h2 * (c + m) * u1 + a2 * e1 * g * m *
(-h1 * l + I1 * u1)) + d2 * (d1 * e1 * I2 * (c + m) * u1 -
(d1 + g) * I2 * (c + m) * u1 + e1 * g * m * (h1 * l - I1 *
u1)) * u2 + a1 * d2 * e1 * l * (c + m) * (h2 * l - I2 * u2)))/(a2 *
e2 * h2 * l * (c + m) * (a1 * e1 * l + (d1 - d1 * e1 + g) *
u1) - a2 * (a1 * e1 * e2 * I2 * l * (c + m) + e2 * (d1 + g) *
I2 * (c + m) * u1 - e1 * (d1 * e2 * I2 * (c + m) * u1 + g *
m * (h1 * l - I1 * u1))) * u2) + (d2 * ((e1 * g * m * (-h1 *
l + I1 * u1) * u2)/((c + m) * (a1 * e1 * l + (d1 - d1 * e1 +
g) * u1)) + e2 * (-h2 * l + I2 * u2)))/(a2 * e2 * l - d2 *
(-1 + e2) * u2), (c * e1 * g * (-h1 * l + I1 * u1))/((c + m) *
(a1 * e1 * l + (d1 - d1 * e1 + g) * u1)), ((d1 + g) * ((d1 -
d1 * e1 + g) * h1 + a1 * e1 * I1) * u1)/(a1 * e1 * (a1 * e1 *
l + (d1 - d1 * e1 + g) * u1)), (e1 * (d1 + g) * (-h1 * l +
I1 * u1))/(a1 * e1 * l + (d1 - d1 * e1 + g) * u1), ((-a1 *
e1 * (d2 * (-1 + e2) * h2 - a2 * e2 * I2) * l * (c + m) + d2 *
(-1 + e2) * (d1 * (-1 + e1) - g) * h2 * (c + m) * u1 + a2 *
(e2 * (d1 + g) * I2 * (c + m) * u1 - e1 * (d1 * e2 * I2 * (c +
m) * u1 + g * m * (h1 * l - I1 * u1)))) * u2 * (l * (-d2 *
(d1 * (-1 + e1) - g) * h2 * (c + m) * u1 + a2 * e1 * g * m *
(-h1 * l + I1 * u1)) + d2 * (d1 * e1 * I2 * (c + m) * u1 -
(d1 + g) * I2 * (c + m) * u1 + e1 * g * m * (h1 * l - I1 *
u1)) * u2 + a1 * d2 * e1 * l * (c + m) * (h2 * l - I2 * u2)))/((c +
m) * (a1 * e1 * l + (d1 - d1 * e1 + g) * u1) * (a2 * e2 * l -
d2 * (-1 + e2) * u2) * (a2 * e2 * h2 * l * (c + m) * (a1 *
e1 * l + (d1 - d1 * e1 + g) * u1) - a2 * (a1 * e1 * e2 * I2 *
l * (c + m) + e2 * (d1 + g) * I2 * (c + m) * u1 - e1 * (d1 *
e2 * I2 * (c + m) * u1 + g * m * (h1 * l - I1 * u1))) * u2)),
(e2 * l * (-a1 * d2 * e1 * h2 * l * (c + m) + d2 * (d1 * (-1 +
e1) - g) * h2 * (c + m) * u1 + a2 * e1 * g * m * (h1 * l -
I1 * u1)) + d2 * e2 * (a1 * e1 * I2 * l * (c + m) + (d1 + g) *
I2 * (c + m) * u1 - e1 * (d1 * I2 * (c + m) * u1 + g * m *
(h1 * l - I1 * u1))) * u2)/((c + m) * (a1 * e1 * l + (d1 -
d1 * e1 + g) * u1) * (a2 * e2 * l - d2 * (-1 + e2) * u2)))
}
# calculate meta-ecosystem biomass stocks for each trophic
# compartment, a.k.a., meta-ecosystem productivity
PREDI2$Ntot <- PREDI2$N1 + PREDI2$N2
PREDI2$Ptot <- PREDI2$P1 + PREDI2$P2
PREDI2$Ctot <- PREDI2$C1 + PREDI2$C2
# now calculate recycling flux for the local ecosystem
PREDI2$FLUX_Eco1_check <- PREDI2$FLUX_P1 + PREDI2$FLUX_C1
PREDI2$FLUX_Eco2_check <- PREDI2$FLUX_P2 + PREDI2$FLUX_C2
# and the meta-ecosystem recycling flux
PREDI2$FLUX_Ptot <- PREDI2$FLUX_P1 + PREDI2$FLUX_P2
PREDI2$FLUX_Ctot <- PREDI2$FLUX_C1 + PREDI2$FLUX_C2
PREDI2$MetaEcoFlux <- PREDI2$FLUX_Eco_1 + PREDI2$FLUX_Eco_2 - PREDI2$FLUX_Q
# Finally, calculate the production for each compartment at the
# meta-ecosystem scale
PREDI2 <- dplyr::mutate(PREDI2, PROD_Ptot = PROD_P1 + PROD_P2, .after = "PROD_C2")
PREDI2 <- dplyr::mutate(PREDI2, PROD_Ctot = PROD_C1 + PROD_C2, .after = "PROD_Ptot")
# Equilibrium stocks less than or equal to 0 make no biological
# sense, so we will flag them in the PREDI2 dataset.
PREDI2 <- PREDI2 %>%
dplyr::mutate(., biosense = ifelse(N1 > 0 & P1 > 0 & C1 > 0 & N2 >
0 & P2 > 0 & C2 > 0 & Q > 0, "yes", "no"), .after = MetaEcoFlux) %>%
dplyr::mutate_at(vars(biosense), factor)
Here we perform stability analyses for the model in this scenario where
I1 >> I2. The rationale and code are the same as above.
# create an empty dataframe
MathStabI2 <- NULL
for (i in 1:nrow(PREDI2)) {
# fetch parameter values from the equilibrium simulations df
I1 = PREDI2$I1[i]
I2 = PREDI2$I2[i]
l = PREDI2$l[i]
u1 = PREDI2$u1[i]
u2 = PREDI2$u2[i]
a1 = PREDI2$a1[i]
a2 = PREDI2$a2[i]
h1 = PREDI2$h1[i]
h2 = PREDI2$h2[i]
d1 = PREDI2$d1[i]
d2 = PREDI2$d2[i]
g = PREDI2$g[i]
m = PREDI2$m[i]
c = PREDI2$c[i]
e1 = PREDI2$e1[i]
e2 = PREDI2$e2[i]
# Put in the solutions to the equilibrium values here:
dN1 = PREDI2$N1[i]
dP1 = PREDI2$P1[i]
dC1 = PREDI2$C1[i]
dN2 = PREDI2$N2[i]
dP2 = PREDI2$P2[i]
dC2 = PREDI2$C2[i]
# create the Jacobian from Mathematica NOTE: in transposing from
# Mathematica, the Jacobian is here assembled by row
Jacob <- rbind(c(-l - dP1 * u1, h1 - dN1 * u1, d1, 0, 0, 0), c(dP1 *
u1, -a1 * dC1 - h1 + dN1 * u1, -a1 * dP1, 0, 0, 0), c(0, a1 * dC1 *
e1, -d1 - g + a1 * e1 * dP1, 0, 0, 0), c(0, 0, 0, -l - dP2 * u2,
h2 - dN2 * u2, d2), c(0, 0, 0, dP2 * u2, -a2 * dC2 - h2 + dN2 *
u2, -a2 * dP2), c(0, 0, (g * m)/(c + m), 0, a2 * dC2 * e2, -d2 +
a2 * e2 * dP2))
# Check to see if the real part of the leading eigenvalue is less
# than 0 or not - if it is than the system is stable.
if (max(Re(base::eigen(Jacob)$values)) < 0) {
stable <- "stable"
} else {
stable <- "unstable"
}
MathStabI2 <- rbind(MathStabI2, data.frame(TIME = i, I1, I2, l, u1,
u2, a1, a2, h1, h2, d1, d2, g, m, c, e1, e2, dN1, dP1, dC1, dN2,
dP2, dC2, EiV1 = Re(eigen(Jacob)$values[1]), EiV2 = Re(eigen(Jacob)$values[2]),
EiV3 = Re(eigen(Jacob)$values[3]), EiV4 = Re(eigen(Jacob)$values[4]),
EiV5 = Re(eigen(Jacob)$values[5]), EiV6 = Re(eigen(Jacob)$values[6]),
maxEv = max(Re(base::eigen(Jacob)$values)), stable = stable, biosense = PREDI2$biosense[i]))
}
MathStabI2$stable <- as.factor(MathStabI2$stable)
# separate unstable equilibria to work with later
MathStabI2US <- subset(MathStabI2, MathStabI2$stable == "unstable")
3.24% of the stability analyses runs produced unstable results (i.e., 324 out of 10000 iterations). Now we check if these are the same ones that produce equilibria where state variables are < 0 at equilibrium (i.e., not biologicaly plausible).
I2biononsense <- subset(PREDI2, PREDI2$biosense == "no")
identical(as.numeric(MathStabI2US[, "TIME"]), I2biononsense[, "TIME"])
[1] TRUE
Since it appears that they are, let’s add the stable/unstable information to
PREDI2.
Now, let’s exclude the unstable parameter sets from the analyses below and from
our results. Then, we will sample 1000 iterations out of PREDI2 to use in
later figures.
PREDI2pos <- filter(PREDI2, PREDI2$biosense == "yes" & PREDI2$stable ==
"stable")
# sample PREDpos to only use 1000 random simulation results
PREDI2_sample1000 <- PREDI2pos[sample(nrow(PREDI2pos), size = 1000), ]
PREDI2_1k <- droplevels(PREDI2_sample1000)
# pivot the dataframe
PREDI2_biomass_long <- pivot_longer(PREDI2_1k, names_to = "Compartment",
values_to = "Stock", cols = c(N1, P1, C1, N2, P2, C2, Q, Ntot, Ptot,
Ctot)) %>%
dplyr::mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
"P1" | Compartment == "C1", "Donor", if_else(Compartment == "N2" |
Compartment == "P2" | Compartment == "C2", "Recipient", if_else(Compartment ==
"Q", "Dispersers' Pool", "Meta-ecosystem")))) %>%
dplyr::mutate(., Ecosystem = fct_relevel(Ecosystem, "Donor", "Recipient",
"Meta-ecosystem")) %>%
dplyr::filter(., Ecosystem != "Dispersers' Pool")
PREDI2_biomass_long$Compartment <- as_factor(PREDI2_biomass_long$Compartment)
PREDI2_biomass_long$Ecosystem <- as_factor(PREDI2_biomass_long$Ecosystem)
comparts_medians <- PREDI2_biomass_long %>%
dplyr::group_by(Compartment, Ecosystem) %>%
dplyr::summarise(., median = median(Stock), .groups = "keep")
I2bm <- ggplot(data = PREDI2_biomass_long, aes(x = Compartment, y = Stock,
col = Ecosystem)) + geom_boxjitter(alpha = 0.25, outlier.intersect = TRUE,
outlier.shape = 25, jitter.shape = 21) + stat_summary(geom = "label",
fun = quantile, fun.args = list(probs = 0.5), aes(label = round(after_stat(y),
2)), position = position_dodge(width = 0.2), show.legend = F, size = 3,
label.padding = unit(0.15, "lines"), label.size = 0.15, fontface = "bold",
color = "black") + facet_grid(. ~ Ecosystem, scales = "free") + scale_color_manual(values = met.brewer("Egypt",
4)) + ggtitle("Organic Biomass and Nutrients Stock") + theme_pubr() +
coord_cartesian(y = c(0, 100)) + theme(legend.position = "none")
I2bm
Figure 6: Effects of along-gradient, diffusive consumers movement in a meta-ecosystem. Compared to the control scenario, the release of primary producers in ecosystem 1 is larger following movement of consumers out of this ecosystem (confront this figure with Figure 1). Environmental fertility values: I1 = 18, I2 = 2. All other specifications as in Figure 1.
FluxPREDI2 <- PREDI2_1k %>%
dplyr::select(., TIME:c, FLUX_P1:FLUX_Eco_2, MetaEcoFlux) %>%
tidyr::pivot_longer(., names_to = "Compartment", values_to = "Flux",
cols = c(FLUX_P1:MetaEcoFlux)) %>%
dplyr::mutate(., Scale = if_else(Compartment == "FLUX_P1" | Compartment ==
"FLUX_C1" | Compartment == "FLUX_Eco_1", "Donor", if_else(Compartment ==
"FLUX_P2" | Compartment == "FLUX_C2" | Compartment == "FLUX_Eco_2",
"Recipient", "Meta-ecosystem"))) %>%
dplyr::mutate(., Scale = fct_relevel(Scale, "Donor", "Recipient", "Meta-ecosystem"))
FluxPREDI2$Scale <- as.factor(FluxPREDI2$Scale)
FluxPREDI2$Compartment <- as.factor(FluxPREDI2$Compartment)
comparts_medians <- FluxPREDI2 %>%
dplyr::group_by(Compartment, Scale) %>%
dplyr::summarise(., median = median(Flux), .groups = "keep")
I2nf <- FluxPREDI2 %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, c("FLUX_P1",
"FLUX_C1", "FLUX_Eco_1", "FLUX_P2", "FLUX_C2", "FLUX_Eco_2", "MetaEcoFlux"))) %>%
ggplot(aes(x = Compartment, y = Flux)) + geom_boxjitter(outlier.intersect = TRUE,
alpha = 0.1, outlier.shape = 25, jitter.shape = 21) + stat_summary(geom = "label",
fun = quantile, fun.args = list(probs = 0.5), aes(label = round(after_stat(y),
2)), position = position_dodge(width = 0.2), show.legend = F, size = 3,
label.padding = unit(0.15, "lines"), label.size = 0.15, fontface = "bold",
color = "black") + facet_grid(. ~ Scale, scales = "free") + coord_cartesian(y = c(0,
300)) + scale_x_discrete(labels = c(FLUX_C1 = "C1", FLUX_P1 = "P1",
FLUX_Eco_1 = "Donor", FLUX_C2 = "C2", FLUX_P2 = "P2", FLUX_Eco_2 = "Recipient",
MetaEcoFlux = "Meta-ecosystem")) + ggtitle("Nutrient Flux") + theme_pubr() +
theme(legend.position = "none", axis.title.x = element_blank())
I2nf
Figure 7: With a higher fertility in the ecosystem 1 (i.e., I1 >> I2) nutrient flux in this ecosystem increases, most likely as a consequence of the release primary producers. Conversely, ecosystem 2 shows lower and less variable levels of nutrient flux. At the meta-ecosystem scale, nutrient flux is higher than in either local ecosystems but lower than in the case of equal environmental fertility across ecosystems (Figure 5). Environmental fertility values: I1 = 18, I2 = 2. All other specifications as in Figure 5.
# I2all <- I2bm + I2nf + plot_annotation(tag_levels = 'a', tag_prefix
# = '(', tag_suffix = ')') + plot_layout(ncol = 1, nrow = 2)
# ggsave(I2all, filename = '../Results/LowI2BmNf.pdf', device =
# 'pdf', dpi = 300, width = 8, height = 6)
The code chunk below analyses the experimental scenario of higher environmental fertility in ecosystem 2. In this case, consumer movement happens counter to the resource availability gradient (i.e., non-diffusive movement; Gounand et al. (2018)).
# Here we set values for inorganic Nutrients input rate in each patch
# In this case, I1 < I2
I1 = 2
I2 = 18
# leaching rate
l = 0.1
# Allocate an empty data frame to save the data in
PREDI1 <- NULL
# Allocate some rows and columns. We will replace the current numbers
# with results below Data we are calculating are (bio)mass for soil
# nutrients (N1, N2), producers (P1, P2), consumers (C1,C2),
# disperses (Q)
PREDI1 <- rbind(PREDI1, data.frame(TIME = seq(0, 100, 1), I1 = I1, I2 = I2,
l = l, u1 = seq(0, 100, 1), u2 = seq(0, 100, 1), a1 = seq(0, 100, 1),
a2 = seq(0, 100, 1), h1 = seq(0, 100, 1), h2 = seq(0, 100, 1), d1 = seq(0,
100, 1), d2 = seq(0, 100, 1), g = seq(0, 100, 1), m = seq(0, 100,
1), c = seq(0, 100, 1), e1 = seq(0, 100, 1), e2 = seq(0, 100, 1),
N1 = seq(0, 100, 1), P1 = seq(0, 100, 1), C1 = seq(0, 100, 1), N2 = seq(0,
100, 1), P2 = seq(0, 100, 1), C2 = seq(0, 100, 1), Q = seq(0, 100,
1), FLUX_P1 = seq(0, 100, 1), FLUX_C1 = seq(0, 100, 1), FLUX_P2 = seq(0,
100, 1), FLUX_C2 = seq(0, 100, 1), FLUX_Eco_1 = seq(0, 100, 1),
FLUX_Eco_2 = seq(0, 100, 1), FLUX_Q = seq(0, 100, 1), PROD_P1 = seq(0,
100, 1), PROD_C1 = seq(0, 100, 1), PROD_P2 = seq(0, 100, 1), PROD_C2 = seq(0,
100, 1)))
# We want to start simulations at 0
i1 = 0
for (i in seq(0, 9999, 1)) {
# go to the next row
i1 = i1 + 1
# Let's sample the parameter values from DATA. We will sample
# these x times as defined by the i for loop above producers
# uptake rates
u1 = DATA1[i1, 2] # in ecosystem 1
u2 = DATA1[i1, 3] # in ecosystem 2
# consumers attack rates
a1 = DATA1[i1, 4]
a2 = DATA1[i1, 5]
# death rates (= recycling rate)
h1 = DATA1[i1, 6] # death (=recycling) rate for producers in ecosystems 1
h2 = DATA1[i1, 7] # death (=recycling) rate for producers in ecosystem 2
d1 = DATA1[i1, 8] # death (=recycling) rate for consumers in ecosystem 1
d2 = DATA1[i1, 9] # death (=recycling) rate for consumers in ecosystem 2
c = DATA1[i1, 10] # death rate in the disperser's pool Q
# movement rate to the Disperser's pool Q
g = DATA1[i1, 11]
# movement rate from the Disperser's pool Q
m = DATA1[i1, 12]
# efficiency of consumers
e1 = DATA2[i1, 2] # in ecosystem 1
e2 = DATA2[i1, 3] # in ecosystem 2
PREDI1[i1, ] <- c(i1, I1, I2, l, u1, u2, a1, a2, h1, h2, d1, d2, g,
m, c, e1, e2, ((d1 - d1 * e1 + g) * h1 + a1 * e1 * I1)/(a1 * e1 *
l + (d1 - d1 * e1 + g) * u1), (d1 + g)/(a1 * e1), (e1 * (-h1 *
l + I1 * u1))/(a1 * e1 * l + (d1 - d1 * e1 + g) * u1), (-a1 *
e1 * (d2 * (-1 + e2) * h2 - a2 * e2 * I2) * l * (c + m) + d2 *
(-1 + e2) * (d1 * (-1 + e1) - g) * h2 * (c + m) * u1 + a2 *
(e2 * (d1 + g) * I2 * (c + m) * u1 - e1 * (d1 * e2 * I2 * (c +
m) * u1 + g * m * (h1 * l - I1 * u1))))/((c + m) * (a1 *
e1 * l + (d1 - d1 * e1 + g) * u1) * (a2 * e2 * l - d2 * (-1 +
e2) * u2)), (l * (-d2 * (d1 * (-1 + e1) - g) * h2 * (c + m) *
u1 + a2 * e1 * g * m * (-h1 * l + I1 * u1)) + d2 * (d1 * e1 *
I2 * (c + m) * u1 - (d1 + g) * I2 * (c + m) * u1 + e1 * g *
m * (h1 * l - I1 * u1)) * u2 + a1 * d2 * e1 * l * (c + m) *
(h2 * l - I2 * u2))/(a2 * e2 * h2 * l * (c + m) * (a1 * e1 *
l + (d1 - d1 * e1 + g) * u1) - a2 * (a1 * e1 * e2 * I2 * l *
(c + m) + e2 * (d1 + g) * I2 * (c + m) * u1 - e1 * (d1 * e2 *
I2 * (c + m) * u1 + g * m * (h1 * l - I1 * u1))) * u2), ((e1 *
g * m * (-h1 * l + I1 * u1) * u2)/((c + m) * (a1 * e1 * l +
(d1 - d1 * e1 + g) * u1)) + e2 * (-h2 * l + I2 * u2))/(a2 *
e2 * l - d2 * (-1 + e2) * u2), (e1 * g * (-h1 * l + I1 * u1))/((c +
m) * (a1 * e1 * l + (d1 - d1 * e1 + g) * u1)), ((d1 + g) *
h1)/(a1 * e1), (d1 * e1 * (-h1 * l + I1 * u1))/(a1 * e1 * l +
(d1 - d1 * e1 + g) * u1), (h2 * (l * (-d2 * (d1 * (-1 + e1) -
g) * h2 * (c + m) * u1 + a2 * e1 * g * m * (-h1 * l + I1 *
u1)) + d2 * (d1 * e1 * I2 * (c + m) * u1 - (d1 + g) * I2 *
(c + m) * u1 + e1 * g * m * (h1 * l - I1 * u1)) * u2 + a1 *
d2 * e1 * l * (c + m) * (h2 * l - I2 * u2)))/(a2 * e2 * h2 *
l * (c + m) * (a1 * e1 * l + (d1 - d1 * e1 + g) * u1) - a2 *
(a1 * e1 * e2 * I2 * l * (c + m) + e2 * (d1 + g) * I2 * (c +
m) * u1 - e1 * (d1 * e2 * I2 * (c + m) * u1 + g * m * (h1 *
l - I1 * u1))) * u2), (d2 * ((e1 * g * m * (-h1 * l + I1 *
u1) * u2)/((c + m) * (a1 * e1 * l + (d1 - d1 * e1 + g) * u1)) +
e2 * (-h2 * l + I2 * u2)))/(a2 * e2 * l - d2 * (-1 + e2) *
u2), ((d1 + g) * h1)/(a1 * e1) + (d1 * e1 * (-h1 * l + I1 *
u1))/(a1 * e1 * l + (d1 - d1 * e1 + g) * u1), (h2 * (l * (-d2 *
(d1 * (-1 + e1) - g) * h2 * (c + m) * u1 + a2 * e1 * g * m *
(-h1 * l + I1 * u1)) + d2 * (d1 * e1 * I2 * (c + m) * u1 -
(d1 + g) * I2 * (c + m) * u1 + e1 * g * m * (h1 * l - I1 *
u1)) * u2 + a1 * d2 * e1 * l * (c + m) * (h2 * l - I2 * u2)))/(a2 *
e2 * h2 * l * (c + m) * (a1 * e1 * l + (d1 - d1 * e1 + g) *
u1) - a2 * (a1 * e1 * e2 * I2 * l * (c + m) + e2 * (d1 + g) *
I2 * (c + m) * u1 - e1 * (d1 * e2 * I2 * (c + m) * u1 + g *
m * (h1 * l - I1 * u1))) * u2) + (d2 * ((e1 * g * m * (-h1 *
l + I1 * u1) * u2)/((c + m) * (a1 * e1 * l + (d1 - d1 * e1 +
g) * u1)) + e2 * (-h2 * l + I2 * u2)))/(a2 * e2 * l - d2 *
(-1 + e2) * u2), (c * e1 * g * (-h1 * l + I1 * u1))/((c + m) *
(a1 * e1 * l + (d1 - d1 * e1 + g) * u1)), ((d1 + g) * ((d1 -
d1 * e1 + g) * h1 + a1 * e1 * I1) * u1)/(a1 * e1 * (a1 * e1 *
l + (d1 - d1 * e1 + g) * u1)), (e1 * (d1 + g) * (-h1 * l +
I1 * u1))/(a1 * e1 * l + (d1 - d1 * e1 + g) * u1), ((-a1 *
e1 * (d2 * (-1 + e2) * h2 - a2 * e2 * I2) * l * (c + m) + d2 *
(-1 + e2) * (d1 * (-1 + e1) - g) * h2 * (c + m) * u1 + a2 *
(e2 * (d1 + g) * I2 * (c + m) * u1 - e1 * (d1 * e2 * I2 * (c +
m) * u1 + g * m * (h1 * l - I1 * u1)))) * u2 * (l * (-d2 *
(d1 * (-1 + e1) - g) * h2 * (c + m) * u1 + a2 * e1 * g * m *
(-h1 * l + I1 * u1)) + d2 * (d1 * e1 * I2 * (c + m) * u1 -
(d1 + g) * I2 * (c + m) * u1 + e1 * g * m * (h1 * l - I1 *
u1)) * u2 + a1 * d2 * e1 * l * (c + m) * (h2 * l - I2 * u2)))/((c +
m) * (a1 * e1 * l + (d1 - d1 * e1 + g) * u1) * (a2 * e2 * l -
d2 * (-1 + e2) * u2) * (a2 * e2 * h2 * l * (c + m) * (a1 *
e1 * l + (d1 - d1 * e1 + g) * u1) - a2 * (a1 * e1 * e2 * I2 *
l * (c + m) + e2 * (d1 + g) * I2 * (c + m) * u1 - e1 * (d1 *
e2 * I2 * (c + m) * u1 + g * m * (h1 * l - I1 * u1))) * u2)),
(e2 * l * (-a1 * d2 * e1 * h2 * l * (c + m) + d2 * (d1 * (-1 +
e1) - g) * h2 * (c + m) * u1 + a2 * e1 * g * m * (h1 * l -
I1 * u1)) + d2 * e2 * (a1 * e1 * I2 * l * (c + m) + (d1 + g) *
I2 * (c + m) * u1 - e1 * (d1 * I2 * (c + m) * u1 + g * m *
(h1 * l - I1 * u1))) * u2)/((c + m) * (a1 * e1 * l + (d1 -
d1 * e1 + g) * u1) * (a2 * e2 * l - d2 * (-1 + e2) * u2)))
}
# calculate meta-ecosystem biomass stocks for each trophic
# compartment, a.k.a., meta-ecosystem productivity
PREDI1$Ntot <- PREDI1$N1 + PREDI1$N2
PREDI1$Ptot <- PREDI1$P1 + PREDI1$P2
PREDI1$Ctot <- PREDI1$C1 + PREDI1$C2
# now calculate recycling flux for the local ecosystem
PREDI1$FLUX_Eco1_check <- PREDI1$FLUX_P1 + PREDI1$FLUX_C1
PREDI1$FLUX_Eco2_check <- PREDI1$FLUX_P2 + PREDI1$FLUX_C2
# and the meta-ecosystem recycling flux
PREDI1$FLUX_Ptot <- PREDI1$FLUX_P1 + PREDI1$FLUX_P2
PREDI1$FLUX_Ctot <- PREDI1$FLUX_C1 + PREDI1$FLUX_C2
PREDI1$MetaEcoFlux <- PREDI1$FLUX_Eco_1 + PREDI1$FLUX_Eco_2 - PREDI1$FLUX_Q
# Finally, calculate the compartment production at meta-ecosystem
# scale
PREDI1 <- dplyr::mutate(PREDI1, PROD_Ptot = PROD_P1 + PROD_P2, .after = "PROD_C2")
PREDI1 <- dplyr::mutate(PREDI1, PROD_Ctot = PROD_C1 + PROD_C2, .after = "PROD_Ptot")
# Equilibrium stocks less than or equal to 0 make no biological
# sense, so we will remove them from the PRED dataset.
PREDI1 <- PREDI1 %>%
dplyr::mutate(., biosense = ifelse(N1 > 0 & P1 > 0 & C1 > 0 & N2 >
0 & P2 > 0 & C2 > 0 & Q > 0, "yes", "no"), .after = MetaEcoFlux) %>%
dplyr::mutate_at(vars(biosense), factor)
Here we perform stability analyses for the scenario of higher environmental
fertility in ecosystem 2, the recipient (i.e., I1 << I2). The analyses
follow the established workflow (see above).
# create an empty dataframe
MathStabI1 <- NULL
for (i in 1:nrow(PREDI1)) {
# fetch parameter values from the equilibrium simulations df
I1 = PREDI1$I1[i]
I2 = PREDI1$I2[i]
l = PREDI1$l[i]
u1 = PREDI1$u1[i]
u2 = PREDI1$u2[i]
a1 = PREDI1$a1[i]
a2 = PREDI1$a2[i]
h1 = PREDI1$h1[i]
h2 = PREDI1$h2[i]
d1 = PREDI1$d1[i]
d2 = PREDI1$d2[i]
g = PREDI1$g[i]
m = PREDI1$m[i]
c = PREDI1$c[i]
e1 = PREDI1$e1[i]
e2 = PREDI1$e2[i]
# Put in the solutions to the equilibrium values here:
dN1 = PREDI1$N1[i]
dP1 = PREDI1$P1[i]
dC1 = PREDI1$C1[i]
dN2 = PREDI1$N2[i]
dP2 = PREDI1$P2[i]
dC2 = PREDI1$C2[i]
# create the Jacobian from Mathematica NOTE: in transposing from
# Mathematica, the Jacobian is here assembled by row
Jacob <- rbind(c(-l - dP1 * u1, h1 - dN1 * u1, d1, 0, 0, 0), c(dP1 *
u1, -a1 * dC1 - h1 + dN1 * u1, -a1 * dP1, 0, 0, 0), c(0, a1 * dC1 *
e1, -d1 - g + a1 * e1 * dP1, 0, 0, 0), c(0, 0, 0, -l - dP2 * u2,
h2 - dN2 * u2, d2), c(0, 0, 0, dP2 * u2, -a2 * dC2 - h2 + dN2 *
u2, -a2 * dP2), c(0, 0, (g * m)/(c + m), 0, a2 * dC2 * e2, -d2 +
a2 * e2 * dP2))
# Check to see if the real part of the leading eigenvalue is less
# than 0 or not - if it is than the system is stable.
if (max(Re(base::eigen(Jacob)$values)) < 0) {
stable <- "stable"
} else {
stable <- "unstable"
}
MathStabI1 <- rbind(MathStabI1, data.frame(TIME = i, I1, I2, l, u1,
u2, a1, a2, h1, h2, d1, d2, g, m, c, e1, e2, dN1, dP1, dC1, dN2,
dP2, dC2, EiV1 = Re(eigen(Jacob)$values[1]), EiV2 = Re(eigen(Jacob)$values[2]),
EiV3 = Re(eigen(Jacob)$values[3]), EiV4 = Re(eigen(Jacob)$values[4]),
EiV5 = Re(eigen(Jacob)$values[5]), EiV6 = Re(eigen(Jacob)$values[6]),
maxEv = max(Re(base::eigen(Jacob)$values)), stable = stable, biosense = PREDI1$biosense[i]))
}
MathStabI1$stable <- as.factor(MathStabI1$stable)
# separate unstable equilibria to work with later
MathStabI1US <- subset(MathStabI1, MathStabI1$stable == "unstable")
2.85% of the stability analyses runs produced unstable results (i.e., 285 out of 10000 iterations). As always, let’s check whether these are also the parameter sets that produce state variable values at equilibrium that do not have biological meaning (i.e., <0).
I1biononsense <- subset(PREDI1, PREDI1$biosense == "no")
identical(as.numeric(MathStabI1US[, "TIME"]), I1biononsense[, "TIME"])
[1] TRUE
Now, since that proved to be the case, let’s add the stable/unstable information to
PREDI1.
Now, let’s exclude the unstable, biological nonsense parameter sets from the
analyses below and from our results. Then, we sample 1000 random iterations
from PREDI1 and store them for later use in figures and comparisons.
PREDI1pos <- filter(PREDI1, PREDI1$biosense == "yes" & PREDI1$stable ==
"stable")
# sample PREDpos to only use 1000 random simulation results
PREDI1_sample1000 <- PREDI1pos[sample(nrow(PREDI1pos), size = 1000), ]
PREDI1_1k <- droplevels(PREDI1_sample1000)
# pivot the
PREDI1_biomass_long <- pivot_longer(PREDI1_1k, names_to = "Compartment",
values_to = "Stock", cols = c(N1, P1, C1, N2, P2, C2, Q, Ntot, Ptot,
Ctot)) %>%
dplyr::mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
"C1" | Compartment == "P1", "Donor", if_else(Compartment == "N2" |
Compartment == "P2" | Compartment == "C2", "Recipient", if_else(Compartment ==
"Q", "Dispersers' Pool", "Meta-ecosystem")))) %>%
dplyr::mutate(., Ecosystem = fct_relevel(Ecosystem, "Donor", "Recipient",
"Meta-ecosystem")) %>%
dplyr::filter(., Ecosystem != "Dispersers' Pool")
PREDI1_biomass_long$Compartment <- as_factor(PREDI1_biomass_long$Compartment)
PREDI1_biomass_long$Ecosystem <- as_factor(PREDI1_biomass_long$Ecosystem)
comparts_medians <- PREDI1_biomass_long %>%
dplyr::group_by(Compartment, Ecosystem) %>%
dplyr::summarise(., median = median(Stock), .groups = "keep")
I1bm <- ggplot(data = PREDI1_biomass_long, aes(x = Compartment, y = Stock,
col = Ecosystem)) + geom_boxjitter(alpha = 0.25, outlier.intersect = TRUE,
outlier.shape = 25, jitter.shape = 21) + stat_summary(geom = "label",
fun = quantile, fun.args = list(probs = 0.5), aes(label = round(after_stat(y),
2)), position = position_dodge(width = 0.2), show.legend = F, size = 3,
label.padding = unit(0.15, "lines"), label.size = 0.15, fontface = "bold",
color = "black") + facet_grid(. ~ Ecosystem, scales = "free") + scale_color_manual(values = met.brewer("Egypt",
4)) + ggtitle("Organic Biomass and Nutrient Stock") + theme_pubr() +
coord_cartesian(y = c(0, 100)) + theme(legend.position = "none")
I1bm
Figure 8: Movement of consumer against the environmental fertility gradient leads to a strong reduction of the nutrients stock in the donor ecosystem 1, stemming from the release of local primary producers from consumer control. In the recipient ecosystem 2, the nutrients stock increases beyond the level of enrichment seen in Figure 1, as local and immigrant consumer strongly suppress the local primary producers. All other specifications as in Figure 1. Environmental fertility values: I1 = 2, I2 = 18.
Following from the stronger trophic cascades, nutrient flux is also affected at both local and meta-ecosystem scales when consumer movement happens against the fertility gradient. Here, median nutrient flux in the recipient, and more fertile, ecosystem 2 is more than double the median flux in the case of diffusive consumer movement (Figure 7). Furthermore, it is also higher than in the case of equal local environmental fertility—i.e., lack of a fertility gradient (Figure 5). Likewise, median nutrient flux at the meta-ecosystem scale is higher than in either the along-gradient (Figure 7) and no-gradient (Figure 5) cases.
FluxPREDI1 <- PREDI1_1k %>%
dplyr::select(TIME:c, FLUX_P1:FLUX_Eco_2, MetaEcoFlux) %>%
tidyr::pivot_longer(names_to = "Compartment", values_to = "Flux", cols = c(FLUX_P1:MetaEcoFlux)) %>%
dplyr::mutate(Scale = if_else(Compartment == "FLUX_P1" | Compartment ==
"FLUX_C1" | Compartment == "FLUX_Eco_1", "Donor", if_else(Compartment ==
"FLUX_P2" | Compartment == "FLUX_C2" | Compartment == "FLUX_Eco_2",
"Recipient", "Meta-ecosystem"))) %>%
dplyr::mutate(., Scale = fct_relevel(Scale, "Donor", "Recipient", "Meta-ecosystem"))
FluxPREDI1$Scale <- as.factor(FluxPREDI1$Scale)
FluxPREDI1$Compartment <- as.factor(FluxPREDI1$Compartment)
comparts_medians <- FluxPREDI1 %>%
dplyr::group_by(Compartment, Scale) %>%
dplyr::summarise(., median = median(Flux), .groups = "keep")
I1nf <- FluxPREDI1 %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, c("FLUX_P1",
"FLUX_C1", "FLUX_Eco_1", "FLUX_P2", "FLUX_C2", "FLUX_Eco_2", "MetaEcoFlux"))) %>%
ggplot(., aes(x = Compartment, y = Flux)) + geom_boxjitter(outlier.intersect = TRUE,
alpha = 0.1, outlier.shape = 25, jitter.shape = 21) + stat_summary(geom = "label",
fun = quantile, fun.args = list(probs = 0.5), aes(label = round(after_stat(y),
2)), position = position_dodge(width = 0.2), show.legend = F, size = 3,
label.padding = unit(0.15, "lines"), label.size = 0.15, fontface = "bold",
color = "black") + facet_grid(. ~ Scale, scales = "free") + theme_pubr() +
coord_cartesian(y = c(0, 400)) + scale_x_discrete(labels = c(FLUX_C1 = "C1",
FLUX_P1 = "P1", FLUX_Eco_1 = "Donor", FLUX_C2 = "C2", FLUX_P2 = "P2",
FLUX_Eco_2 = "Recipient", MetaEcoFlux = "Meta-ecosystem")) + ggtitle("Nutrient Flux") +
theme_pubr() + theme(legend.position = "none", axis.title.x = element_blank())
I1nf
Figure 9: Nutrient flux more than doubles in ecosystem 2, when this is the more fertile of the two ecosystem and is also the recipient of consumer movement. In turn, this leads to an overall increase in the nutrient flux at the meta-ecosystem scale, compared to scenarios where environmental fertility is equal among local ecosystem (Figure 5) or is higher in the donor ecosystem (Figure 19). Environmental fertility values: I1 = 2, I2 = 18.
# I1all <- I1bm + I1nf + plot_annotation(tag_levels = 'a', tag_prefix
# = '(', tag_suffix = ')') + plot_layout(ncol = 1, nrow = 2)
# ggsave(I1all, filename = '../Results/LowI1BmNf.pdf', device =
# 'pdf', dpi = 300, width = 8, height = 6)
Here, we summarize the results of the model’s stability analyses, and collect this information to produce Table S1 shown in the Supplementary Materials. As the table shows, in all three scenario tested, only a small portion of the parameter sets used to simulate the model’s behavior at equilibrium produces unstable results. As detailed above, we exclude these unstable parameter sets from further investigations.
The table is saved in folder ../Results/ as a .tex file.
# summarise the number of stable and unstable parameter sets for each
# stability analysis run with the numerical approximation method and
# calculate the percentage of unstable parameter sets over the total
# iterations of the simulations
MStabSum <- MathStab %>%
dplyr::group_by(., stable) %>%
dplyr::summarise(., n = n(), .groups = "keep") %>%
add_column(., Fertility = "Control") %>%
pivot_wider(id_cols = Fertility, names_from = stable, values_from = n) %>%
dplyr::mutate(., Percent = ((unstable/nrow(MathStab)) * 100))
I2MStabSum <- MathStabI2 %>%
dplyr::group_by(., stable) %>%
dplyr::summarise(., n = n(), .groups = "keep") %>%
add_column(., Fertility = "Along-gradient consumer movement") %>%
pivot_wider(id_cols = Fertility, names_from = stable, values_from = n) %>%
dplyr::mutate(., Percent = ((unstable/nrow(MathStabI2)) * 100))
I1MStabSum <- MathStabI1 %>%
dplyr::group_by(., stable) %>%
dplyr::summarise(., n = n(), .groups = "keep") %>%
add_column(., Fertility = "Against-gradient consumer movement") %>%
pivot_wider(id_cols = Fertility, names_from = stable, values_from = n) %>%
dplyr::mutate(., Percent = ((unstable/nrow(MathStabI1)) * 100))
# Bind the objects generated above into a single dataframe and then
# produce a table that shows the number of stable and unstable
# parameter sets for each stability analysis method
SAsummary <- bind_rows(MStabSum, I2MStabSum, I1MStabSum) %>%
ungroup() %>%
gt() %>%
fmt_number(columns = 4, decimals = 2) %>%
cols_label(., Fertility = "Fertility", stable = md("Stable"), unstable = md("Unstable"),
Percent = md("\\%")) %>%
tab_header(title = md("**Summary of stability analyses**, showing the number of
Stable and Unstable equilibria, and the percentage of Unstable
equilibria over the total number of iterations (n = 10000).
Rows are grouped by consumer movement scenario.")) %>%
tab_style(style = list(cell_text(weight = "bold")), locations = cells_column_labels()) %>%
tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())
# save it gtsave(SAsummary, filename =
# '../Results/StabilityAnalyses_SummaryTable.tex')
# view it
SAsummary
| Summary of stability analyses, showing the number of Stable and Unstable equilibria, and the percentage of Unstable equilibria over the total number of iterations (n = 10000). Rows are grouped by consumer movement scenario. | |||
| Fertility | Stable | Unstable | % |
|---|---|---|---|
| Control | 9845 | 155 | 1.55 |
| Along-gradient consumer movement | 9676 | 324 | 3.24 |
| Against-gradient consumer movement | 9715 | 285 | 2.85 |
Here we summarize the above results and produce Figures 2, 3, and 4 shown in the main text and Figures S1 to S5 in the Supplementary Materials. We refer the reader to those documents to see the actual figures; here we only present the code used to generate them.
Note that Figures 2, 3, and 4 in the main text use a log10-transformed response ratio (\(LRR_X\), where \(X \in [N, P, C]\)) to assess change in nutrient stock and biomass accumulation, nutrient flux, and trophic compartment productivity in different consumer movement scenarios. Our reference for calculating the \(LRR_X\) are the results for the control scenario (Control scenario: gradient-neutral consumer movement), where environmental fertility is homogeneous between local ecosystems (i.e., I1 = I2 = 10).
Mathematically, the ratio’s formula is:
\(\displaystyle LRR_X = log_{10} \biggl(\frac{X_{i, I_1 \gtrless I_2}}{X_{i, I_1 = I_2}}\biggr)\)
where \(X \in [N, P, C]\) is the trophic compartment of interest, \(i\) is either the local or meta-ecosystem, and \(I_1 \gtrless I_2\) at the numerator represents the direction of fertility gradient in the experimental scenario of interest. At the denominator, \(I_1 = I_2\) represents the equal fertility conditions of the control scenario. Please, see the Numerical Analyses section of the manuscript as well as Supplementary Figure A.1 for a primer on how to interpret the LRR graphs.
# Here we create a custom ggplot2 theme to use in the graphs below
consmov_theme <- theme_pubr() + theme(legend.position = "bottom", legend.text = element_text(size = 16),
legend.title = element_text(size = 17), plot.title = element_text(size = 17),
axis.title.x = element_text(size = 17, face = "bold"), axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 15), axis.title.y = element_text(size = 17,
face = "bold"), strip.text = element_text(size = 16))
We are working with the PRED, PREDI2, and PREDI1 datasets we produced
above (sections Control scenario: gradient-neutral consumer movement and
Experimental scenarios: along- and against-gradient consumer movement).
First, we will check whether there is overlap among the equilibria that are
either unstable or without biological meaning across PRED, PREDI2, and
PREDI1. Then, we will join them into a new object, EnvFert_res.
Mode FALSE TRUE
logical 61 94
Mode FALSE TRUE
logical 29 126
There is no overlap among the three datasets in terms of equilibria that are unstable and lack biological meaning. This is something that will need accounting for later on when we work with the LRR. For now, we can remove all equilibria that are either unstable or do not have biological sense as we are going to create a figure using the raw data.
EnvFert_respos <- filter(EnvFert_res, biosense == "yes")
We now transform EnvFert_respos to a long-format dataframe and add two
categorical variables: Function and Scale. Function tells us which
ecosystem function is measured. Scale shows instead whether the function
is being measured at the local or meta-ecosystem scale.
EnvFert_long <- select(EnvFert_respos, !u1:e2) %>%
pivot_longer(., cols = N1:PROD_Ctot, names_to = "EcoFunction", values_to = "Value") %>%
dplyr::mutate(., Function = if_else(EcoFunction == "N1" | EcoFunction ==
"N2" | EcoFunction == "P1" | EcoFunction == "P2" | EcoFunction ==
"C1" | EcoFunction == "C2" | EcoFunction == "Ntot" | EcoFunction ==
"Ptot" | EcoFunction == "Ctot", "Stock", if_else(EcoFunction ==
"FLUX_P1" | EcoFunction == "FLUX_P2" | EcoFunction == "FLUX_C1" |
EcoFunction == "FLUX_C2" | EcoFunction == "FLUX_Ptot" | EcoFunction ==
"FLUX_Ctot", "Nutrient Flux", "Trophic Productivity"))) %>%
dplyr::mutate(., Scale = if_else(EcoFunction == "N1" | EcoFunction ==
"P1" | EcoFunction == "C1" | EcoFunction == "FLUX_P1" | EcoFunction ==
"FLUX_C1" | EcoFunction == "PROD_P1" | EcoFunction == "PROD_C1",
"Ecosystem 1", if_else(EcoFunction == "N2" | EcoFunction == "P2" |
EcoFunction == "C2" | EcoFunction == "FLUX_P2" | EcoFunction ==
"FLUX_C2" | EcoFunction == "PROD_P2" | EcoFunction == "PROD_C2",
"Ecosystem 2", "Meta-ecosystem"))) %>%
dplyr::mutate(., Fertility = factor(Fertility), Function = factor(Function),
Scale = factor(Scale), EcoFunction = factor(EcoFunction))
head(EnvFert_long)
# A tibble: 6 × 9
biosense Fertility TIME I1 I2 EcoFunction Value Funct…¹ Scale
<fct> <fct> <dbl> <dbl> <dbl> <fct> <dbl> <fct> <fct>
1 yes Equal 1 10 10 N1 3.07 Stock Ecos…
2 yes Equal 1 10 10 P1 1.43 Stock Ecos…
3 yes Equal 1 10 10 C1 1.31 Stock Ecos…
4 yes Equal 1 10 10 N2 3.26 Stock Ecos…
5 yes Equal 1 10 10 P2 1.55 Stock Ecos…
6 yes Equal 1 10 10 C2 2.68 Stock Ecos…
# … with abbreviated variable name ¹Function
Now that we have a dataset ready for plotting, we produce the graphs. The next code chunk produces Figure A.2 in the Supplementary Materials. Note that we show the code, but not the figure itself as it is presented in the Supplementary Materials.
StockSumm <- EnvFert_long %>%
filter(., Function == "Stock") %>%
dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, Nutrients = "N1",
Nutrients = "N2", Nutrients = "Ntot", `Primary\n Producers` = "P1",
`Primary\n Producers` = "P2", `Primary\n Producers` = "Ptot", Consumers = "C1",
Consumers = "C2", Consumers = "Ctot")) %>%
dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Nutrients",
"Primary\n Producers", "Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along",
"Against")) %>%
ggplot(., aes(x = EcoFunction, y = Value)) + geom_hline(yintercept = 10,
linetype = "dashed", lwd = 0.25) + geom_boxplot(aes(color = Fertility)) +
scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Gradient-neutral",
"Along-gradient", "Against-gradient")) + facet_grid(. ~ Scale,
scales = "free") + labs(title = "Nutrient Stock and Biomass") + ylab(" ") +
xlab(" ") + coord_cartesian(ylim = c(0, 75))
FluxSumm <- EnvFert_long %>%
filter(., Function == "Nutrient Flux") %>%
dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, `Primary\n Producers` = "FLUX_P1",
`Primary\n Producers` = "FLUX_P2", `Primary\n Producers` = "FLUX_Ptot",
Consumers = "FLUX_C1", Consumers = "FLUX_C2", Consumers = "FLUX_Ctot")) %>%
dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along",
"Against")) %>%
ggplot(., aes(x = EcoFunction, y = Value)) + geom_hline(yintercept = 50,
linetype = "dashed", lwd = 0.25) + geom_boxplot(aes(color = Fertility)) +
scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Gradient-neutral",
"Along-gradient", "Against-gradient")) + facet_grid(. ~ Scale,
scales = "free") + labs(title = "Nutrient Flux") + ylab(" ") + xlab(" ") +
coord_cartesian(ylim = c(0, 350))
ProdSumm <- EnvFert_long %>%
filter(., Function == "Trophic Productivity") %>%
dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, `Primary\n Producers` = "PROD_P1",
`Primary\n Producers` = "PROD_P2", `Primary\n Producers` = "PROD_Ptot",
Consumers = "PROD_C1", Consumers = "PROD_C2", Consumers = "PROD_Ctot")) %>%
dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along",
"Against")) %>%
ggplot(., aes(x = EcoFunction, y = Value)) + geom_hline(yintercept = 50,
linetype = "dashed", lwd = 0.25) + geom_boxplot(aes(color = Fertility)) +
scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Gradient-neutral",
"Along-gradient", "Against-gradient")) + facet_grid(. ~ Scale,
scales = "free") + labs(title = "Trophic Compartment Productivity") +
ylab(" ") + xlab("Trophic compartment") + coord_cartesian(ylim = c(0,
500))
FuncAll <- StockSumm + FluxSumm + ProdSumm + plot_layout(ncol = 1, nrow = 3,
guides = "collect") + plot_annotation(tag_levels = "a", tag_prefix = "(",
tag_suffix = ")") & consmov_theme
# FuncAll
# ggsave(FuncAll, filename = '../Results/FunctionSummary_10kN.pdf',
# device = 'pdf', dpi = 300, width = 11, height = 10)
Now, we produce Figures 2, 3, and 4 from the main text. First, we are going
to create new dataframes (PREDI2rr and PREDI1rr1) that contain only the
response ratio values for nutrient stock and biomass accumulation, nutrient flux,
and trophic compartment productivity of all state variables. These dataframes
will also contain information about the fertility conditions under which we
measure each ecosystem function. After calculating the response ratio, we will
remove the unstable, biologically not meaningful equilibria from
PREDI1 and PREDI2 before merging them in a new object, EnvFertRR.
In calculating the LRR, the difference in the number of equilibria that are
either unstable or lack biological meaning among PRED, PREDI1, and PREDI2
matters. To account for this, let us first create two objects, I1NsUsrm and
I2NsUsrm that contain the unstable, nonsensical equilibria that are found in
PRED but are not included in either PREDI1 and PREDI2. We are going to
remove the ones in PREDI1 and PREDI2 making use of the columns biosense
and stable. Conversely, we will use I1NsUsrm and I2NsUsrm to remove
instance of the response ratio being calculated using a denominator that is
from an unstable, nonsensical equilibrium.
Then, let’s create two new objects, PREDI2rr and PREDI1rr, to store the
response ratio values calculated by dividing ecosystem function values in
PREDI2 and PREDI1 by those in PRED, respectively. We will then merge these
two objects in a single one called EnvFertRR, and calculate the \(C{:}R\) biomass
ratio.
# need to keep raw data in the same df as the ratio to select the >0
# cases for state variables
PREDI2rr <- PREDI2 %>%
# we are going to batch-create new response ratio versions of
# each variable we start with the 'stock' group of variables
dplyr::mutate(., N1rr = N1/PRED$N1, P1rr = P1/PRED$P1, C1rr = C1/PRED$C1,
N2rr = N2/PRED$N2, P2rr = P2/PRED$P2, C2rr = C2/PRED$C2, Nrr = Ntot/PRED$Ntot,
Prr = Ptot/PRED$Ptot, Crr = Ctot/PRED$Ctot) %>%
# next, we work on the 'flux' group of variables
dplyr::mutate(., FLUX_P1rr = FLUX_P1/PRED$FLUX_P1, FLUX_C1rr = FLUX_C1/PRED$FLUX_C1,
FLUX_P2rr = FLUX_P2/PRED$FLUX_P2, FLUX_C2rr = FLUX_C2/PRED$FLUX_C2,
FLUX_Prr = FLUX_Ptot/PRED$FLUX_Ptot, FLUX_Crr = FLUX_Ctot/PRED$FLUX_Ctot) %>%
# finally, we work on the 'productivity' group of variables
dplyr::mutate(., PROD_P1rr = PROD_P1/PRED$PROD_P1, PROD_C1rr = PROD_C1/PRED$PROD_C1,
PROD_P2rr = PROD_P2/PRED$PROD_P2, PROD_C2rr = PROD_C2/PRED$PROD_C2,
PROD_Prr = PROD_Ptot/PRED$PROD_Ptot, PROD_Crr = PROD_Ctot/PRED$PROD_Ctot) %>%
# and now let's subset the dataframe to include the resp ratio
# columns
dplyr::select(., TIME:I2, N1:Ctot, N1rr:PROD_Crr, biosense, stable, u1:e2) %>%
dplyr::filter(., biosense == "yes" | stable == "yes") %>%
dplyr::filter(., !(TIME %in% I2NsUsrm$TIME))
# repeat, but for PREDI1
PREDI1rr <- PREDI1 %>%
# we are going to batch-create new response ratio versions of
# each variable we start with the 'stock' group of variables
dplyr::mutate(., N1rr = N1/PRED$N1, P1rr = P1/PRED$P1, C1rr = C1/PRED$C1,
N2rr = N2/PRED$N2, P2rr = P2/PRED$P2, C2rr = C2/PRED$C2, Nrr = Ntot/PRED$Ntot,
Prr = Ptot/PRED$Ptot, Crr = Ctot/PRED$Ctot) %>%
# next, we work on the 'flux' group of variables
dplyr::mutate(., FLUX_P1rr = FLUX_P1/PRED$FLUX_P1, FLUX_C1rr = FLUX_C1/PRED$FLUX_C1,
FLUX_P2rr = FLUX_P2/PRED$FLUX_P2, FLUX_C2rr = FLUX_C2/PRED$FLUX_C2,
FLUX_Prr = FLUX_Ptot/PRED$FLUX_Ptot, FLUX_Crr = FLUX_Ctot/PRED$FLUX_Ctot) %>%
# finally, we work on the 'productivity' group of variables
dplyr::mutate(., PROD_P1rr = PROD_P1/PRED$PROD_P1, PROD_C1rr = PROD_C1/PRED$PROD_C1,
PROD_P2rr = PROD_P2/PRED$PROD_P2, PROD_C2rr = PROD_C2/PRED$PROD_C2,
PROD_Prr = PROD_Ptot/PRED$PROD_Ptot, PROD_Crr = PROD_Ctot/PRED$PROD_Ctot) %>%
# and now let's subset the dataframe to include only the resp
# ratio columns
dplyr::select(., TIME:I2, N1:Ctot, N1rr:PROD_Crr, biosense, stable, u1:e2) %>%
dplyr::filter(., biosense == "yes" | stable == "yes") %>%
dplyr::filter(., !(TIME %in% I1NsUsrm$TIME))
EnvFertRR <- bind_rows(Along = PREDI2rr, Against = PREDI1rr, .id = "Fertility") %>%
dplyr::select(., Fertility, biosense, stable, TIME:I2, N1rr:PROD_Crr,
N1:Ctot, u1:e2) %>%
dplyr::mutate(., CR_Eco1 = C1/P1, CR_Eco2 = C2/P2, CR_MetaEco = Ctot/Ptot)
# now build the graph, first by pivoting to longer format, then
# graphing
Note that object EnvFertRR contains 19 301 and not 20
000 as would be expected from merging PREDI1 and PREDI2, because of the two
rounds of removal of unstable equilibria after said merge. The first round
removes those unstable, biologically not meaningful equilibria contained in
PREDI1 and PREDI2, the second uses objects I1NsUsrm and I2NsUsrm to
remove the ones that are contained in PRED but do not overlap with those in
PREDI1 or PREDI2.
Now, we build the graphs. First, pivot the EnvFertRR dataframe from
wide to long format.
EnvFertRR_long <- select(EnvFertRR, !(N1:e2)) %>%
pivot_longer(., cols = N1rr:PROD_Crr, names_to = "Compartment", values_to = "Ratio") %>%
dplyr::mutate(., Function = if_else(Compartment == "FLUX_P1rr" | Compartment ==
"FLUX_P2rr" | Compartment == "FLUX_C1rr" | Compartment == "FLUX_C2rr" |
Compartment == "FLUX_Prr" | Compartment == "FLUX_Crr", "Nutrient Flux",
if_else(Compartment == "PROD_P1rr" | Compartment == "PROD_C1rr" |
Compartment == "PROD_P2rr" | Compartment == "PROD_C2rr" | Compartment ==
"PROD_Prr" | Compartment == "PROD_Crr", "Trophic Productivity",
"Stock"))) %>%
dplyr::mutate(., Scale = if_else(Compartment == "N1rr" | Compartment ==
"P1rr" | Compartment == "C1rr" | Compartment == "FLUX_P1rr" | Compartment ==
"FLUX_C1rr" | Compartment == "PROD_P1rr" | Compartment == "PROD_C1rr",
"Ecosystem 1", if_else(Compartment == "N2rr" | Compartment == "P2rr" |
Compartment == "C2rr" | Compartment == "FLUX_P2rr" | Compartment ==
"FLUX_C2rr" | Compartment == "PROD_P2rr" | Compartment == "PROD_C2rr",
"Ecosystem 2", "Meta-ecosystem"))) %>%
dplyr::mutate(., Fertility = factor(Fertility), Function = factor(Function),
Scale = factor(Scale), Compartment = factor(Compartment)) %>%
dplyr::mutate(., biosense = fct_drop(biosense), stable = fct_drop(stable))
In the following code chunks, we subset and then plot each ecosystem function individually against the LRR to create Figure 2, 3, and 4 in the main text. Note, again, that here we show only the code used to generate the figures but not the figures themselves as these can be found in the main text.
The following code chunk produces Figure 2, showing changes in nutrient stocks and biomass accumulation among experimental and control scenarios.
ptHCsubset <- c("#DDAA33", "#BB5566")
ptHCsubset_fill <- c("#f8eed6", "#f1dde0")
EFstockLRR <- EnvFertRR_long %>%
filter(., Function == "Stock") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, Nutrients = "N1rr",
Nutrients = "N2rr", Nutrients = "Nrr", `Primary\n Producers` = "P1rr",
`Primary\n Producers` = "P2rr", `Primary\n Producers` = "Prr",
Consumers = "C1rr", Consumers = "C2rr", Consumers = "Crr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Nutrients",
"Primary\n Producers", "Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
ggplot(., aes(x = Compartment, y = log10(Ratio))) + geom_rect(aes(xmin = -Inf,
xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") + geom_hline(yintercept = 0,
linetype = "dashed", colour = "black") + geom_boxplot(aes(color = Fertility,
fill = Fertility)) + scale_color_manual(values = ptHCsubset, name = "Consumer movement",
labels = c("Along-gradient", "Against-gradient")) + scale_fill_manual(values = ptHCsubset_fill,
name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
facet_grid(. ~ Scale, scales = "free") + labs(title = "Nutrient Stock and Biomass") +
ylab("LRR") + xlab("Trophic compartment") + consmov_theme
# ggsave(EFstockLRR, filename = '../Results/Stock_LRR_10kN.pdf',
# device = 'pdf', dpi = 300, width = 10, height = 6)
The following code chunk produces Figure 3, which shows changes to nutrient flux among experimental and control scenarios.
EFfluxLRR <- EnvFertRR_long %>%
filter(., Function == "Nutrient Flux") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, `Primary\n Producers` = "FLUX_P1rr",
`Primary\n Producers` = "FLUX_P2rr", `Primary\n Producers` = "FLUX_Prr",
Consumers = "FLUX_C1rr", Consumers = "FLUX_C2rr", Consumers = "FLUX_Crr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
ggplot(., aes(x = Compartment, y = log10(Ratio))) + geom_rect(aes(xmin = -Inf,
xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") + geom_hline(yintercept = 0,
linetype = "dashed", colour = "black") + geom_boxplot(aes(color = Fertility,
fill = Fertility)) + scale_color_manual(values = ptHCsubset, name = "Consumer movement",
labels = c("Along-gradient", "Against-gradient")) + scale_fill_manual(values = ptHCsubset_fill,
name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
facet_grid(. ~ Scale, scales = "free") + labs(title = "Nutrient Flux") +
ylab("LRR") + xlab("Trophic compartment") + consmov_theme
# ggsave(EFfluxLRR, filename = '../Results/Flux_LRR_10kN.pdf', device
# = 'pdf', dpi = 300, width = 10, height = 6)
And, finally, here we produce Figure 4, that shows change in the trophic compartment productivity among experimental and control scenarios.
EFprodLRR <- EnvFertRR_long %>%
filter(., Function == "Trophic Productivity") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, `Primary\n Producers` = "PROD_P1rr",
`Primary\n Producers` = "PROD_P2rr", `Primary\n Producers` = "PROD_Prr",
Consumers = "PROD_C1rr", Consumers = "PROD_C2rr", Consumers = "PROD_Crr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
ggplot(., aes(x = Compartment, y = log10(Ratio))) + geom_rect(aes(xmin = -Inf,
xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") + geom_hline(yintercept = 0,
linetype = "dashed", colour = "black") + geom_boxplot(aes(color = Fertility,
fill = Fertility)) + scale_color_manual(values = ptHCsubset, name = "Consumer movement",
labels = c("Along-gradient", "Against-gradient")) + scale_fill_manual(values = ptHCsubset_fill,
name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
facet_grid(. ~ Scale, scales = "free") + labs(title = "Trophic Compartment Productivity") +
ylab("LRR") + xlab("Trophic compartment") + consmov_theme
# ggsave(EFprodLRR, filename = '../Results/Production_LRR_10kN.pdf',
# device = 'pdf', dpi = 300, width = 10, height = 6)
In the two sections below, we investigate if our results are robust to changes in biomass distribution within the system caused by moving consumers and produce Figures S2 to S5 in the Supplementary Materials.
Here, we focus on those parameter sets that produce \(C{:}R < 1\). That is, on those parameter sets that produce classic biomass pyramids (McCauley et al. 2018).
As before, we begin with graphs for the untransformed data and then proceed to the LRR graphs.
First, using the object EnvFert_CRbelow10 created above that exclude all \(C{:}R\)
values above 10, we further exclude all \(C{:}R\) values above 1 and create a new
object, EnvFert_CRbelow1.
Now, we lengthen this new dataset for plotting with ggplot2.
EnvFert_CRbelow1_long <- select(EnvFert_CRbelow1, !u1:e2) %>%
pivot_longer(., cols = N1:CR_MetaEco, names_to = "EcoFunction", values_to = "Value") %>%
dplyr::mutate(., Function = if_else(EcoFunction == "N1" | EcoFunction ==
"N2" | EcoFunction == "P1" | EcoFunction == "P2" | EcoFunction ==
"C1" | EcoFunction == "C2" | EcoFunction == "Ntot" | EcoFunction ==
"Ptot" | EcoFunction == "Ctot", "Stock", if_else(EcoFunction ==
"FLUX_P1" | EcoFunction == "FLUX_P2" | EcoFunction == "FLUX_C1" |
EcoFunction == "FLUX_C2" | EcoFunction == "FLUX_Ptot" | EcoFunction ==
"FLUX_Ctot", "Nutrient Flux", ifelse(EcoFunction == "CR_Eco1" |
EcoFunction == "CR_Eco2" | EcoFunction == "CR_MetaEco", "C:R",
"Trophic Productivity")))) %>%
dplyr::mutate(., Scale = if_else(EcoFunction == "N1" | EcoFunction ==
"P1" | EcoFunction == "C1" | EcoFunction == "FLUX_P1" | EcoFunction ==
"FLUX_C1" | EcoFunction == "PROD_P1" | EcoFunction == "PROD_C1" |
EcoFunction == "CR_Eco1", "Ecosystem 1", if_else(EcoFunction ==
"N2" | EcoFunction == "P2" | EcoFunction == "C2" | EcoFunction ==
"FLUX_P2" | EcoFunction == "FLUX_C2" | EcoFunction == "PROD_P2" |
EcoFunction == "PROD_C2" | EcoFunction == "CR_Eco2", "Ecosystem 2",
"Meta-ecosystem"))) %>%
dplyr::mutate(., Fertility = factor(Fertility), Function = factor(Function),
Scale = factor(Scale), EcoFunction = factor(EcoFunction))
head(EnvFert_CRbelow1_long)
# A tibble: 6 × 9
biosense Fertility TIME I1 I2 EcoFunct…¹ Value Funct…² Scale
<fct> <fct> <dbl> <dbl> <dbl> <fct> <dbl> <fct> <fct>
1 yes Equal 2 10 10 N1 0.773 Stock Ecos…
2 yes Equal 2 10 10 P1 16.7 Stock Ecos…
3 yes Equal 2 10 10 C1 0.429 Stock Ecos…
4 yes Equal 2 10 10 N2 0.858 Stock Ecos…
5 yes Equal 2 10 10 P2 12.4 Stock Ecos…
6 yes Equal 2 10 10 C2 3.20 Stock Ecos…
# … with abbreviated variable names ¹EcoFunction, ²Function
And now we use the newly created dataframe EnvFert_CRbelow1_long to produce a
figure analogous to Supplementary Figure A.2. Unlike previous figures, the
following one does not appear in the Supplementary Material and is thus shown here.
StockSummCRbelow1 <- EnvFert_CRbelow1_long %>%
filter(., Function == "Stock") %>%
dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, Nutrients = "N1",
Nutrients = "N2", Nutrients = "Ntot", `Primary\n Producers` = "P1",
`Primary\n Producers` = "P2", `Primary\n Producers` = "Ptot", Consumers = "C1",
Consumers = "C2", Consumers = "Ctot")) %>%
dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Nutrients",
"Primary\n Producers", "Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along",
"Against")) %>%
ggplot(., aes(x = EcoFunction, y = Value)) + geom_hline(yintercept = 10,
linetype = "dashed", lwd = 0.25) + geom_boxplot(aes(color = Fertility)) +
scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Equal",
"Along-gradient", "Against-gradient")) + facet_grid(. ~ Scale,
scales = "free") + labs(title = "Nutrient Stock and Biomass") + ylab(" ") +
xlab(" ") + coord_cartesian(ylim = c(0, 75))
FluxSummCRbelow1 <- EnvFert_CRbelow1_long %>%
filter(., Function == "Nutrient Flux") %>%
dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, `Primary\n Producers` = "FLUX_P1",
`Primary\n Producers` = "FLUX_P2", `Primary\n Producers` = "FLUX_Ptot",
Consumers = "FLUX_C1", Consumers = "FLUX_C2", Consumers = "FLUX_Ctot")) %>%
dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along",
"Against")) %>%
ggplot(., aes(x = EcoFunction, y = Value)) + geom_hline(yintercept = 50,
linetype = "dashed", lwd = 0.25) + geom_boxplot(aes(color = Fertility)) +
scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Equal",
"Along-gradient", "Against-gradient")) + facet_grid(. ~ Scale,
scales = "free") + labs(title = "Nutrient Flux") + ylab(" ") + xlab(" ") +
coord_cartesian(ylim = c(0, 350))
ProdSummCRbelow1 <- EnvFert_CRbelow1_long %>%
filter(., Function == "Trophic Productivity") %>%
dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, `Primary\n Producers` = "PROD_P1",
`Primary\n Producers` = "PROD_P2", `Primary\n Producers` = "PROD_Ptot",
Consumers = "PROD_C1", Consumers = "PROD_C2", Consumers = "PROD_Ctot")) %>%
dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along",
"Against")) %>%
ggplot(., aes(x = EcoFunction, y = Value)) + geom_hline(yintercept = 50,
linetype = "dashed", lwd = 0.25) + geom_boxplot(aes(color = Fertility)) +
scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Equal",
"Along-gradient", "Against-gradient")) + facet_grid(. ~ Scale,
scales = "free") + labs(title = "Trophic Compartment Productivity") +
ylab(" ") + xlab("Trophic Compartment") + coord_cartesian(ylim = c(0,
500))
FuncAllCRbelow1 <- StockSummCRbelow1 + FluxSummCRbelow1 + ProdSummCRbelow1 +
plot_layout(ncol = 1, nrow = 3, guides = "collect") + plot_annotation(tag_levels = "a",
tag_prefix = "(", tag_suffix = ")") & consmov_theme
FuncAllCRbelow1
Figure 10: Change in nutrient stock and biomass (a), nutrient flux (b), and trophic compartment productivity (c) when consumers move in a heterogeneous meta-ecosystem, considering only parameters sets for which C:R < 1 (unstransformed data). Notably, we observe a clear quantative difference in the nutrient stock at the meta-ecosystem scale, compared with Figure A.2. However, this does not appear to cause a qualitative change in the patterns of local and meta-ecosystem functioning observe, as noted in Figure A.3. All specifications as in Figure A.2.
While the quantitative changes are more noticeable here, again in the distance between the upper and lower hinges of the boxes, there does not seem to be a marked change in the qualitative pattern observed in our results for any of the ecosystem functions of interest.
Let’s proceed to work on the LRR values.
As we just did for the untransformed data, let’s use the previously created
object EnvFertRR to select the parameter sets producing \(C{:}R < 1\)
and store them in a new object, EnvFertRR_CRbelow1.
EnvFertRR_CRbelow1 <- dplyr::filter(EnvFertRR, CR_Eco1 <= 1 & CR_Eco2 <=
1 & CR_MetaEco <= 1)
We now lengthen this newly created object.
EnvFertRR_CRbelow1_long <- select(EnvFertRR_CRbelow1, !(N1:e2)) %>%
pivot_longer(., cols = N1rr:CR_MetaEco, names_to = "Compartment", values_to = "Ratio") %>%
dplyr::mutate(., Function = if_else(Compartment == "FLUX_P1rr" | Compartment ==
"FLUX_P2rr" | Compartment == "FLUX_C1rr" | Compartment == "FLUX_C2rr" |
Compartment == "FLUX_Prr" | Compartment == "FLUX_Crr", "Nutrient Flux",
if_else(Compartment == "PROD_P1rr" | Compartment == "PROD_C1rr" |
Compartment == "PROD_P2rr" | Compartment == "PROD_C2rr" | Compartment ==
"PROD_Prr" | Compartment == "PROD_Crr", "Trophic Productivity",
ifelse(Compartment == "CR_Eco1" | Compartment == "CR_Eco2" |
Compartment == "CR_MetaEco", "C:R", "Stock")))) %>%
dplyr::mutate(., Scale = if_else(Compartment == "N1rr" | Compartment ==
"P1rr" | Compartment == "C1rr" | Compartment == "FLUX_P1rr" | Compartment ==
"FLUX_C1rr" | Compartment == "PROD_P1rr" | Compartment == "PROD_C1rr" |
Compartment == "CR_Eco1", "Ecosystem 1", if_else(Compartment ==
"N2rr" | Compartment == "P2rr" | Compartment == "C2rr" | Compartment ==
"FLUX_P2rr" | Compartment == "FLUX_C2rr" | Compartment == "PROD_P2rr" |
Compartment == "PROD_C2rr" | Compartment == "CR_Eco2", "Ecosystem 2",
"Meta-ecosystem"))) %>%
dplyr::mutate(., Fertility = factor(Fertility), Function = factor(Function),
Scale = factor(Scale), Compartment = factor(Compartment)) %>%
dplyr::mutate(., biosense = fct_drop(biosense), stable = fct_drop(stable))
And finally we produce Supplementary Figure A.3. Accordingly, we show the code here but not the figure itself.
# we make sure the y-axis of this figure has the same range of the
# 'all data' graph above by extracting the limits of the y-axis from
# that graph and using them in this one
LRRstock_ymin <- layer_scales(EFstockLRR)$y$get_limits()[1]
LRRstock_ymax <- layer_scales(EFstockLRR)$y$get_limits()[2]
EFstockLRR_CRbelow1 <- EnvFertRR_CRbelow1_long %>%
filter(., Function == "Stock") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, Nutrients = "N1rr",
Nutrients = "N2rr", Nutrients = "Nrr", `Primary\n Producers` = "P1rr",
`Primary\n Producers` = "P2rr", `Primary\n Producers` = "Prr",
Consumers = "C1rr", Consumers = "C2rr", Consumers = "Crr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Nutrients",
"Primary\n Producers", "Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
ggplot(., aes(x = Compartment, y = log10(Ratio))) + geom_rect(aes(xmin = -Inf,
xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") + geom_hline(yintercept = 0,
linetype = "dashed", colour = "black") + geom_boxplot(aes(color = Fertility,
fill = Fertility)) + scale_color_manual(values = ptHCsubset, name = "Consumer movement",
labels = c("Along-gradient", "Against-gradient")) + scale_fill_manual(values = ptHCsubset_fill,
name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
facet_grid(. ~ Scale, scales = "free") + labs(title = "Nutrient Stock and Biomass") +
ylab("LRR") + xlab(" ") + coord_cartesian(ylim = c(LRRstock_ymin, LRRstock_ymax))
# again, we grab the y-axis limits from the relevant all-data graph
LRRflux_ymin <- layer_scales(EFfluxLRR)$y$get_limits()[1]
LRRflux_ymax <- layer_scales(EFfluxLRR)$y$get_limits()[2]
EFfluxLRR_CRbelow1 <- EnvFertRR_CRbelow1_long %>%
filter(., Function == "Nutrient Flux") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, `Primary\n Producers` = "FLUX_P1rr",
`Primary\n Producers` = "FLUX_P2rr", `Primary\n Producers` = "FLUX_Prr",
Consumers = "FLUX_C1rr", Consumers = "FLUX_C2rr", Consumers = "FLUX_Crr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
ggplot(., aes(x = Compartment, y = log10(Ratio))) + geom_rect(aes(xmin = -Inf,
xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") + geom_hline(yintercept = 0,
linetype = "dashed", colour = "black") + geom_boxplot(aes(color = Fertility,
fill = Fertility)) + scale_color_manual(values = ptHCsubset, name = "Consumer movement",
labels = c("Along-gradient", "Against-gradient")) + scale_fill_manual(values = ptHCsubset_fill,
name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
facet_grid(. ~ Scale, scales = "free") + labs(title = "Nutrient Flux") +
ylab("LRR") + xlab(" ") + coord_cartesian(ylim = c(LRRflux_ymin, LRRflux_ymax))
# again, we grab the y-axis limits from the relevant all-data graph
LRRprod_ymin <- layer_scales(EFprodLRR)$y$get_limits()[1]
LRRprod_ymax <- layer_scales(EFprodLRR)$y$get_limits()[2]
EFprodLRR_CRbelow1 <- EnvFertRR_CRbelow1_long %>%
filter(., Function == "Trophic Productivity") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, `Primary\n Producers` = "PROD_P1rr",
`Primary\n Producers` = "PROD_P2rr", `Primary\n Producers` = "PROD_Prr",
Consumers = "PROD_C1rr", Consumers = "PROD_C2rr", Consumers = "PROD_Crr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
ggplot(., aes(x = Compartment, y = log10(Ratio))) + geom_rect(aes(xmin = -Inf,
xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") + geom_hline(yintercept = 0,
linetype = "dashed", colour = "black") + geom_boxplot(aes(color = Fertility,
fill = Fertility)) + scale_color_manual(values = ptHCsubset, name = "Consumer movement",
labels = c("Along-gradient", "Against-gradient")) + scale_fill_manual(values = ptHCsubset_fill,
name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
facet_grid(. ~ Scale, scales = "free") + labs(title = "Trophic Compartment Productivity") +
ylab("LRR") + xlab("Trophic compartment") + coord_cartesian(ylim = c(LRRprod_ymin,
LRRprod_ymax))
EFallLRR_CRbelow1 <- EFstockLRR_CRbelow1 + EFfluxLRR_CRbelow1 + EFprodLRR_CRbelow1 +
plot_layout(ncol = 1, nrow = 3, guides = "collect") + plot_annotation(tag_levels = "a",
tag_prefix = "(", tag_suffix = ")") & consmov_theme
# EFallLRR_CRbelow1
# ggsave(EFallLRR_CRbelow1, filename =
# '../Results/FunctionSummary_LRR_CRbelow1.pdf', device = 'pdf', dpi
# = 300, width = 10, height = 10)
And here, we produce Supplementary Figure A.5, showing median parameter values for the subset of results where \(C{:}R < 1\). Again, we show the code but not the figure itself.
# table of median parameter values EnvFertRR_CRbelow1 %>% select(.,
# Fertility, u1:e2) %>% ungroup() %>% group_by(., Fertility) %>%
# pivot_longer(., cols = -Fertility, names_to = 'Parameter',
# values_to = 'value') %>% mutate(., Parameter =
# as_factor(Parameter)) %>% ungroup() %>% group_by(., Fertility,
# Parameter) %>% summarise(., across(.cols = value, .fns = list(mean
# = mean, sd = sd)), .groups = 'keep') %>% pivot_wider(., names_from
# = 'Fertility', values_from = c('value_mean', 'value_sd'))%>%
# gt(rowname_col = 'Parameter', groupname_col = 'Fertility') %>%
# fmt_number(., columns = 2:5, decimals = 2) %>% tab_spanner(label =
# md('Along'), columns = c('value_mean_Along', 'value_sd_Along')) %>%
# tab_spanner(label = md('Against'), columns =
# c('value_mean_Against', 'value_sd_Against')) %>%
# cols_label(value_mean_Against = 'mean', value_sd_Against = 'sd',
# value_mean_Along = 'mean', value_sd_Along = 'sd')
# %>% gtsave(., filename = '../Results/MeanParam_LRR_CRbelow1.rtf')
EnvFertRR_CRbelow1 %>%
select(., Fertility, u1:e2) %>%
ungroup() %>%
group_by(., Fertility) %>%
pivot_longer(., cols = -Fertility, names_to = "Parameter", values_to = "value") %>%
mutate(., Parameter = as_factor(Parameter)) %>%
ggplot(., aes(x = Parameter, y = value)) + geom_boxjitter(aes(fill = Parameter,
col = Parameter), jitter.fill = NULL, alpha = 0.1) + scale_color_manual(values = met.brewer("Signac",
13)) + scale_fill_manual(values = met.brewer("Signac", 13)) + theme_pubr(legend = "none") +
ylab("Mean") + facet_grid(. ~ Fertility, labeller = as_labeller(c(Against = "Against-gradient Consumer movement",
Along = "Along-gradient Consumer movement")))
# ggsave(filename = '../Results/MeanParam_LRR_CRbelow1.pdf', dpi =
# 300, device = 'pdf')
Here we reproduce Figures 2, 3, and 4 from the main text and Supplementary Figure A.2 using only parameter sets that produce \(C{:}R \in (1,10)\). We exclude \(C{:}R > 10\) values, as these are likely extremely rare in real-world ecosystems. We begin, as usual, with the untransformed dataframe.
First, let’s exclude from EnvFert_respos all the parameter sets producing
\(C{:}R > 10\) and those producing \(C{:}R < 1\).
Now we lengthen the dataset EnvFert_CRabove1 created in the code chunk above
for plotting.
EnvFert_CRabove1_long <- select(EnvFert_CRabove1, !u1:e2) %>%
pivot_longer(., cols = N1:CR_MetaEco, names_to = "EcoFunction", values_to = "Value") %>%
dplyr::mutate(., Function = if_else(EcoFunction == "N1" | EcoFunction ==
"N2" | EcoFunction == "P1" | EcoFunction == "P2" | EcoFunction ==
"C1" | EcoFunction == "C2" | EcoFunction == "Ntot" | EcoFunction ==
"Ptot" | EcoFunction == "Ctot", "Stock", if_else(EcoFunction ==
"FLUX_P1" | EcoFunction == "FLUX_P2" | EcoFunction == "FLUX_C1" |
EcoFunction == "FLUX_C2" | EcoFunction == "FLUX_Ptot" | EcoFunction ==
"FLUX_Ctot", "Nutrient Flux", ifelse(EcoFunction == "CR_Eco1" |
EcoFunction == "CR_Eco2" | EcoFunction == "CR_MetaEco", "C:R",
"Trophic Productivity")))) %>%
dplyr::mutate(., Scale = if_else(EcoFunction == "N1" | EcoFunction ==
"P1" | EcoFunction == "C1" | EcoFunction == "FLUX_P1" | EcoFunction ==
"FLUX_C1" | EcoFunction == "PROD_P1" | EcoFunction == "PROD_C1" |
EcoFunction == "CR_Eco1", "Ecosystem 1", if_else(EcoFunction ==
"N2" | EcoFunction == "P2" | EcoFunction == "C2" | EcoFunction ==
"FLUX_P2" | EcoFunction == "FLUX_C2" | EcoFunction == "PROD_P2" |
EcoFunction == "PROD_C2" | EcoFunction == "CR_Eco2", "Ecosystem 2",
"Meta-ecosystem"))) %>%
dplyr::mutate(., Fertility = factor(Fertility), Function = factor(Function),
Scale = factor(Scale), EcoFunction = factor(EcoFunction))
head(EnvFert_CRabove1_long)
# A tibble: 6 × 9
biosense Fertility TIME I1 I2 EcoFunct…¹ Value Funct…² Scale
<fct> <fct> <dbl> <dbl> <dbl> <fct> <dbl> <fct> <fct>
1 yes Equal 13 10 10 N1 7.29 Stock Ecos…
2 yes Equal 13 10 10 P1 0.848 Stock Ecos…
3 yes Equal 13 10 10 C1 3.26 Stock Ecos…
4 yes Equal 13 10 10 N2 15.5 Stock Ecos…
5 yes Equal 13 10 10 P2 0.894 Stock Ecos…
6 yes Equal 13 10 10 C2 7.23 Stock Ecos…
# … with abbreviated variable names ¹EcoFunction, ²Function
And, finally, we reproduce Supplementary Figure A.2. Note that, again, this figure does not appear in the Supplementary Materials and is thus shown here.
StockSummCRabove1 <- EnvFert_CRabove1_long %>%
filter(., Function == "Stock") %>%
dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, Nutrients = "N1",
Nutrients = "N2", Nutrients = "Ntot", `Primary\n Producers` = "P1",
`Primary\n Producers` = "P2", `Primary\n Producers` = "Ptot", Consumers = "C1",
Consumers = "C2", Consumers = "Ctot")) %>%
dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Nutrients",
"Primary\n Producers", "Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along",
"Against")) %>%
ggplot(., aes(x = EcoFunction, y = Value)) + geom_hline(yintercept = 10,
linetype = "dashed", lwd = 0.25) + geom_boxplot(aes(color = Fertility)) +
scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Equal",
"Along-gradient", "Against-gradient")) + facet_grid(. ~ Scale,
scales = "free") + labs(title = "Nutrient Stock and Biomass") + ylab(" ") +
xlab(" ") + coord_cartesian(ylim = c(0, 75))
FluxSummCRabove1 <- EnvFert_CRabove1_long %>%
filter(., Function == "Nutrient Flux") %>%
dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, `Primary\n Producers` = "FLUX_P1",
`Primary\n Producers` = "FLUX_P2", `Primary\n Producers` = "FLUX_Ptot",
Consumers = "FLUX_C1", Consumers = "FLUX_C2", Consumers = "FLUX_Ctot")) %>%
dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along",
"Against")) %>%
ggplot(., aes(x = EcoFunction, y = Value)) + geom_hline(yintercept = 50,
linetype = "dashed", lwd = 0.25) + geom_boxplot(aes(color = Fertility)) +
scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Equal",
"Along-gradient", "Against-gradient")) + facet_grid(. ~ Scale,
scales = "free") + labs(title = "Nutrient Flux") + ylab(" ") + xlab(" ") +
coord_cartesian(ylim = c(0, 350))
ProdSummCRabove1 <- EnvFert_CRabove1_long %>%
filter(., Function == "Trophic Productivity") %>%
dplyr::mutate(., EcoFunction = fct_recode(EcoFunction, `Primary\n Producers` = "PROD_P1",
`Primary\n Producers` = "PROD_P2", `Primary\n Producers` = "PROD_Ptot",
Consumers = "PROD_C1", Consumers = "PROD_C2", Consumers = "PROD_Ctot")) %>%
dplyr::mutate(., EcoFunction = fct_relevel(EcoFunction, "Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Equal", "Along",
"Against")) %>%
ggplot(., aes(x = EcoFunction, y = Value)) + geom_hline(yintercept = 50,
linetype = "dashed", lwd = 0.25) + geom_boxplot(aes(color = Fertility)) +
scale_color_highcontrast(reverse = F, name = "Consumer movement", labels = c("Equal",
"Along-gradient", "Against-gradient")) + facet_grid(. ~ Scale,
scales = "free") + labs(title = "Trophic Compartment Productivity") +
ylab(" ") + xlab("Trophic compartment") + coord_cartesian(ylim = c(0,
500))
FuncAllCRabove1 <- StockSummCRabove1 + FluxSummCRabove1 + ProdSummCRabove1 +
plot_layout(ncol = 1, nrow = 3, guides = "collect") + plot_annotation(tag_levels = "a",
tag_prefix = "(", tag_suffix = ")") & consmov_theme
FuncAllCRabove1
Figure 11: Change in nutrient stock and biomass (a), nutrient flux (b), and trophic compartment productivity (c) when consumers move in a heterogeneous meta-ecosystem, considering only parameter sets for which C:R > 1 and < 10 (unstransformed data). Here, biomass and nutrient stock results appear to more clearly point towards a trophic cascade at play at the meta-ecosystem scale, compared to Figure A.2. Patterns of nutrient flux and trophic compartment productivity appear qualitatively unchanged. All specifications as in Figure A.2.
# ggsave(FuncAllCRabove1, filename =
# '../Results/FunctionSummary_CRabove1.pdf', device = 'pdf', dpi =
# 300, width = 10, height = 10)
Compared with Supplementary Figure A.2, we see a clear change in the pattern of biomass accumulation—consistent with the inverted, top-heavy biomass pyramid expected when \(C{:}R > 1\) (McCauley et al. 2018).
Here, we produce Supplementary Figure A.4. We use only parameter sets that produce \(C{:}R \in (1,10)\). First, we exclude all \(C{:}R\) values below 1.
Then, we lengthen the dataset for plotting.
EnvFertRR_CRabove1_long <- select(EnvFertRR_CRabove1, !(N1:e2)) %>%
pivot_longer(., cols = N1rr:CR_MetaEco, names_to = "Compartment", values_to = "Ratio") %>%
dplyr::mutate(., Function = if_else(Compartment == "FLUX_P1rr" | Compartment ==
"FLUX_P2rr" | Compartment == "FLUX_C1rr" | Compartment == "FLUX_C2rr" |
Compartment == "FLUX_Prr" | Compartment == "FLUX_Crr", "Nutrient Flux",
if_else(Compartment == "PROD_P1rr" | Compartment == "PROD_C1rr" |
Compartment == "PROD_P2rr" | Compartment == "PROD_C2rr" | Compartment ==
"PROD_Prr" | Compartment == "PROD_Crr", "Trophic Productivity",
ifelse(Compartment == "CR_Eco1" | Compartment == "CR_Eco2" |
Compartment == "CR_MetaEco", "C:R", "Stock")))) %>%
dplyr::mutate(., Scale = if_else(Compartment == "N1rr" | Compartment ==
"P1rr" | Compartment == "C1rr" | Compartment == "FLUX_P1rr" | Compartment ==
"FLUX_C1rr" | Compartment == "PROD_P1rr" | Compartment == "PROD_C1rr" |
Compartment == "CR_Eco1", "Ecosystem 1", if_else(Compartment ==
"N2rr" | Compartment == "P2rr" | Compartment == "C2rr" | Compartment ==
"FLUX_P2rr" | Compartment == "FLUX_C2rr" | Compartment == "PROD_P2rr" |
Compartment == "PROD_C2rr" | Compartment == "CR_Eco2", "Ecosystem 2",
"Meta-ecosystem"))) %>%
dplyr::mutate(., Fertility = factor(Fertility), Function = factor(Function),
Scale = factor(Scale), Compartment = factor(Compartment)) %>%
dplyr::mutate(., biosense = fct_drop(biosense), stable = fct_drop(stable))
And now we produce Supplementary Figure A.4. Note that, as above, we do not show the figure itself here but only the code needed to produce it.
EFstockLRR_CRabove1 <- EnvFertRR_CRabove1_long %>%
filter(., Function == "Stock") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, Nutrients = "N1rr",
Nutrients = "N2rr", Nutrients = "Nrr", `Primary\n Producers` = "P1rr",
`Primary\n Producers` = "P2rr", `Primary\n Producers` = "Prr",
Consumers = "C1rr", Consumers = "C2rr", Consumers = "Crr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Nutrients",
"Primary\n Producers", "Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
ggplot(., aes(x = Compartment, y = log10(Ratio))) + geom_rect(aes(xmin = -Inf,
xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") + geom_hline(yintercept = 0,
linetype = "dashed", colour = "black") + geom_boxplot(aes(color = Fertility,
fill = Fertility)) + scale_color_manual(values = ptHCsubset, name = "Consumer movement",
labels = c("Along-gradient", "Against-gradient")) + scale_fill_manual(values = ptHCsubset_fill,
name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
facet_grid(. ~ Scale, scales = "free") + labs(title = "Nutrient Stock and Biomass") +
ylab("LRR") + xlab(" ") + coord_cartesian(ylim = c(LRRstock_ymin, LRRstock_ymax))
EFfluxLRR_CRabove1 <- EnvFertRR_CRabove1_long %>%
filter(., Function == "Nutrient Flux") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, `Primary\n Producers` = "FLUX_P1rr",
`Primary\n Producers` = "FLUX_P2rr", `Primary\n Producers` = "FLUX_Prr",
Consumers = "FLUX_C1rr", Consumers = "FLUX_C2rr", Consumers = "FLUX_Crr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
ggplot(., aes(x = Compartment, y = log10(Ratio))) + geom_rect(aes(xmin = -Inf,
xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") + geom_hline(yintercept = 0,
linetype = "dashed", colour = "black") + geom_boxplot(aes(color = Fertility,
fill = Fertility)) + scale_color_manual(values = ptHCsubset, name = "Consumer movement",
labels = c("Along-gradient", "Against-gradient")) + scale_fill_manual(values = ptHCsubset_fill,
name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
facet_grid(. ~ Scale, scales = "free") + labs(title = "Nutrient Flux") +
ylab("LRR") + xlab(" ") + coord_cartesian(ylim = c(LRRflux_ymin, LRRflux_ymax))
EFprodLRR_CRabove1 <- EnvFertRR_CRabove1_long %>%
filter(., Function == "Trophic Productivity") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, `Primary\n Producers` = "PROD_P1rr",
`Primary\n Producers` = "PROD_P2rr", `Primary\n Producers` = "PROD_Prr",
Consumers = "PROD_C1rr", Consumers = "PROD_C2rr", Consumers = "PROD_Crr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Primary\n Producers",
"Consumers")) %>%
dplyr::mutate(., Fertility = fct_relevel(Fertility, "Along", "Against")) %>%
ggplot(., aes(x = Compartment, y = log10(Ratio))) + geom_rect(aes(xmin = -Inf,
xmax = Inf, ymin = -Inf, ymax = 0), fill = "gray95") + geom_hline(yintercept = 0,
linetype = "dashed", colour = "black") + geom_boxplot(aes(color = Fertility,
fill = Fertility)) + scale_color_manual(values = ptHCsubset, name = "Consumer movement",
labels = c("Along-gradient", "Against-gradient")) + scale_fill_manual(values = ptHCsubset_fill,
name = "Consumer movement", labels = c("Along-gradient", "Against-gradient")) +
facet_grid(. ~ Scale, scales = "free") + labs(title = "Trophic Compartment Productivity") +
ylab("LRR") + xlab("Trophic compartment") + coord_cartesian(ylim = c(LRRprod_ymin,
LRRprod_ymax))
EFallLRR_CRabove1 <- EFstockLRR_CRabove1 + EFfluxLRR_CRabove1 + EFprodLRR_CRabove1 +
plot_layout(ncol = 1, nrow = 3, guides = "collect") + plot_annotation(tag_levels = "a",
tag_prefix = "(", tag_suffix = ")") & consmov_theme
# ggsave(EFallLRR_CRabove1, filename =
# '../Results/FunctionSummary_LRR_CRabove1.pdf', device = 'pdf', dpi
# = 300, width = 10, height = 10)
As above, we now plot the median parameter values for these parameter sets, which are shown in Supplementary Figure A.6. Accordingly, we only show the code to produce the figure and not the figure itself.
# table of median parameter sets
# EnvFertRR_CRabove1 %>% select(., Fertility, u1:e2) %>% ungroup()
# %>% group_by(., Fertility) %>% pivot_longer(., cols = -Fertility,
# names_to = 'Parameter', values_to = 'value') %>% mutate(.,
# Parameter = as_factor(Parameter)) %>% ungroup() %>% group_by(.,
# Fertility, Parameter) %>% summarise(., across(.cols = value, .fns =
# list(mean = mean, sd = sd)), .groups = 'keep') %>% pivot_wider(.,
# names_from = 'Fertility', values_from = c('value_mean',
# 'value_sd'))%>% gt(rowname_col = 'Parameter', groupname_col =
# 'Fertility') %>% fmt_number(., columns = 2:5, decimals = 2) %>%
# tab_spanner(label = md('Along'), columns = c('value_mean_Along',
# 'value_sd_Along')) %>% tab_spanner(label = md('Against'), columns =
# c('value_mean_Against', 'value_sd_Against')) %>%
# cols_label(value_mean_Against = 'mean', value_sd_Against = 'sd',
# value_mean_Along = 'mean', value_sd_Along = 'sd')
EnvFertRR_CRabove1 %>%
select(., Fertility, u1:e2) %>%
ungroup() %>%
group_by(., Fertility) %>%
pivot_longer(., cols = -Fertility, names_to = "Parameter", values_to = "value") %>%
mutate(., Parameter = as_factor(Parameter)) %>%
ggplot(., aes(x = Parameter, y = value)) + geom_boxjitter(aes(fill = Parameter,
col = Parameter), jitter.fill = NULL, alpha = 0.1) + scale_color_manual(values = met.brewer("Signac",
13)) + scale_fill_manual(values = met.brewer("Signac", 13)) + theme_pubr(legend = "none") +
ylab("Mean") + facet_grid(. ~ Fertility, labeller = as_labeller(c(Against = "Against-gradient Consumer movement",
Along = "Along-gradient Consumer movement")))
# ggsave(filename = '../Results/MeanParam_LRR_CRabove1.pdf', dpi =
# 300, device = 'pdf')
Here, we produce summary tables for the results of the model’s iterations, using the LRR and \(C{:}R\) values we calculated to produce the summary plots above. For LRR, we report median values for each ecosystem function of interest in both local ecosystems and in the meta-ecosystem.
For ease of interpretation, median \(LRR < 0\) will be shown in shades of magenta based on how close they are to -1 and median \(LRR > 0\) will be shown in shades of green based on their closeness to 1.
First, we produce dataframes to contain the median values, one for each of stock, nutrient flux, and primary and secondary productivity.
# first, compute median of log10 ratio values
EnvFertLRR_tables <- dplyr::filter(EnvFertRR_long, Ratio >= 0) %>%
dplyr::group_by(., Fertility, Function, Scale, Compartment) %>%
dplyr::mutate(., LRR = log10(Ratio)) %>%
dplyr::summarise(., min = min(LRR), median = median(LRR), max = max(LRR),
.groups = "keep")
Now, we produce the tables, separating between along and against-gradient movement, as above.
EFStockTab <- EnvFertLRR_tables %>%
dplyr::ungroup() %>%
dplyr::filter(., Function == "Stock") %>%
pivot_wider(., names_from = Fertility, values_from = c(min, median,
max)) %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, Consumers = "C1rr",
Consumers = "C2rr", Consumers = "Crr", Producers = "P1rr", Producers = "P2rr",
Producers = "Prr", Nutrients = "N1rr", Nutrients = "N2rr", Nutrients = "Nrr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Nutrients",
"Producers", "Consumers")) %>%
gt(groupname_col = "Scale") %>%
cols_hide("Function") %>%
tab_spanner(label = md("**Along-gradient**"), columns = c("min_Along",
"median_Along", "max_Along")) %>%
tab_spanner(label = md("**Against-gradient**"), columns = c("min_Against",
"median_Against", "max_Against")) %>%
cols_label(Compartment = md("**Compartment**"), min_Along = md("*Min.*"),
median_Along = md("*Median*"), max_Along = md("*Max.*"), min_Against = md("*Min.*"),
median_Against = md("*Median*"), max_Against = md("*Max.*")) %>%
fmt_number(., columns = c(4:9), decimals = 2) %>%
data_color(., columns = c(6, 7), colors = col_numeric(palette = "PiYG",
domain = c(-1, 1))) %>%
tab_header(title = md("**Change in nutrient stock and organic biomass accumulation** in the meta-ecosystem when consumers move against or along an environmental fertility gradient."),
subtitle = md("Values are expressed as median _LRR_ for 10000 iterations of the model.")) %>%
tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())
# gtsave(EFStockTab, filename = 'EFbiostock.tex', path =
# '../Results/')
EFStockTab
| Change in nutrient stock and organic biomass accumulation in the meta-ecosystem when consumers move against or along an environmental fertility gradient. | ||||||
| Values are expressed as median LRR for 10000 iterations of the model. | ||||||
| Compartment | Against-gradient | Along-gradient | ||||
|---|---|---|---|---|---|---|
| Min. | Median | Max. | Min. | Median | Max. | |
| Ecosystem 1 | ||||||
| Consumers | −3.16 | −0.72 | −0.70 | 0.26 | 0.26 | 2.45 |
| Nutrients | −0.70 | −0.16 | 0.00 | 0.00 | 0.12 | 0.26 |
| Producers | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |
| Ecosystem 2 | ||||||
| Consumers | −0.80 | 0.17 | 0.92 | −1.35 | −0.27 | 0.31 |
| Nutrients | −0.68 | 0.11 | 0.26 | −0.70 | −0.14 | 0.25 |
| Producers | 0.00 | 0.04 | 1.46 | −2.72 | −0.12 | 0.00 |
| Meta-ecosystem | ||||||
| Consumers | −1.16 | 0.08 | 0.59 | −1.08 | −0.09 | 0.57 |
| Nutrients | −0.67 | 0.03 | 0.25 | −0.68 | −0.03 | 0.25 |
| Producers | 0.00 | 0.01 | 0.90 | −0.80 | −0.02 | 0.00 |
EFFluxTab <- EnvFertLRR_tables %>%
dplyr::ungroup() %>%
dplyr::filter(., Function == "Nutrient Flux") %>%
pivot_wider(., names_from = Fertility, values_from = c(min, median,
max)) %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, Consumers = "FLUX_C1rr",
Consumers = "FLUX_C2rr", Consumers = "FLUX_Crr", Producers = "FLUX_P1rr",
Producers = "FLUX_P2rr", Producers = "FLUX_Prr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Producers",
"Consumers")) %>%
gt(groupname_col = "Scale") %>%
cols_hide("Function") %>%
tab_spanner(label = md("**Along-gradient**"), columns = c("min_Along",
"median_Along", "max_Along")) %>%
tab_spanner(label = md("**Against-gradient**"), columns = c("min_Against",
"median_Against", "max_Against")) %>%
cols_label(Compartment = md("**Compartment**"), min_Along = md("*Min.*"),
median_Along = md("*Median*"), max_Along = md("*Max.*"), min_Against = md("*Min.*"),
median_Against = md("*Median*"), max_Against = md("*Max.*")) %>%
fmt_number(., columns = c(4:9), decimals = 2) %>%
data_color(., columns = c(6, 7), colors = col_numeric(palette = "PiYG",
domain = c(-1, 1))) %>%
tab_header(title = md("**Change in nutrient flux** in the meta-ecosystem when consumers move against or along an environmental fertility gradient."),
subtitle = md("Values are expressed as median _LRR_ for 10000 iterations of the model.")) %>%
tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())
# gtsave(EFFluxTab, filename = 'EFnutflux.tex', path = '../Results/')
EFFluxTab
| Change in nutrient flux in the meta-ecosystem when consumers move against or along an environmental fertility gradient. | ||||||
| Values are expressed as median LRR for 10000 iterations of the model. | ||||||
| Compartment | Against-gradient | Along-gradient | ||||
|---|---|---|---|---|---|---|
| Min. | Median | Max. | Min. | Median | Max. | |
| Ecosystem 1 | ||||||
| Consumers | −3.16 | −0.72 | −0.70 | 0.26 | 0.26 | 2.45 |
| Producers | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |
| Ecosystem 2 | ||||||
| Consumers | −0.80 | 0.17 | 0.92 | −1.35 | −0.27 | 0.31 |
| Producers | 0.00 | 0.04 | 1.46 | −2.72 | −0.12 | 0.00 |
| Meta-ecosystem | ||||||
| Consumers | −0.86 | 0.07 | 0.49 | −1.06 | −0.09 | 0.57 |
| Producers | 0.00 | 0.01 | 0.90 | −0.86 | −0.02 | 0.00 |
EFProdTab <- EnvFertLRR_tables %>%
dplyr::ungroup() %>%
dplyr::filter(., Function == "Trophic Productivity") %>%
pivot_wider(., names_from = Fertility, values_from = c(min, median,
max)) %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, Consumers = "PROD_C1rr",
Consumers = "PROD_C2rr", Consumers = "PROD_Crr", Producers = "PROD_P1rr",
Producers = "PROD_P2rr", Producers = "PROD_Prr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Producers",
"Consumers")) %>%
gt(groupname_col = "Scale") %>%
cols_hide("Function") %>%
tab_spanner(label = md("**Along-gradient**"), columns = c("min_Along",
"median_Along", "max_Along")) %>%
tab_spanner(label = md("**Against-gradient**"), columns = c("min_Against",
"median_Against", "max_Against")) %>%
cols_label(Compartment = md("**Compartment**"), min_Along = md("*Min.*"),
median_Along = md("*Median*"), max_Along = md("*Max.*"), min_Against = md("*Min.*"),
median_Against = md("*Median*"), max_Against = md("*Max.*")) %>%
fmt_number(., columns = c(4:9), decimals = 2) %>%
data_color(., columns = c(6, 7), colors = col_numeric(palette = "PiYG",
domain = c(-1, 1))) %>%
tab_header(title = md("**Change in trophic compartment productivity** in the meta-ecosystem when consumers move _along-_ or _against-gradient_ of nutrient availability."),
subtitle = md("Values are expressed as median _LRR_ for 10 000 iterations of the model.")) %>%
tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())
# gtsave(EFProdTab, filename = 'EFprod.tex', path = '../Results/')
EFProdTab
| Change in trophic compartment productivity in the meta-ecosystem when consumers move along- or against-gradient of nutrient availability. | ||||||
| Values are expressed as median LRR for 10 000 iterations of the model. | ||||||
| Compartment | Against-gradient | Along-gradient | ||||
|---|---|---|---|---|---|---|
| Min. | Median | Max. | Min. | Median | Max. | |
| Ecosystem 1 | ||||||
| Consumers | −3.16 | −0.72 | −0.70 | 0.26 | 0.26 | 2.45 |
| Producers | −0.70 | −0.16 | 0.00 | 0.00 | 0.12 | 0.26 |
| Ecosystem 2 | ||||||
| Consumers | 0.01 | 0.22 | 1.40 | −2.63 | −0.47 | −0.01 |
| Producers | 0.00 | 0.18 | 1.40 | −2.63 | −0.34 | 0.00 |
| Meta-ecosystem | ||||||
| Consumers | −0.88 | 0.03 | 0.47 | −1.02 | −0.04 | 0.70 |
| Producers | −0.47 | 0.03 | 0.88 | −0.81 | −0.04 | 0.22 |
Here, we investigate whether different types of consumer movement between the two local ecosystems lead to changes in the distribution of biomass in the meta-ecosystem. Biomass distribution in local and meta-ecosystems can vary from bottom-heavy to top-heavy (McCauley et al. 2018). Bottom-heavy ecosystems present the classic biomass pyramid, with a large base of primary producers supporting an increasingly smaller biomass of primary and secondary consumers. This bottom-heavy distribution is typical of terrestrial ecosystems—albeit exceptions exist (Hatton et al. 2015; McCauley et al. 2018). Conversely, in top-heavy ecosystems, a small base of primary producers supports a large biomass of primary and secondary consumers. This inverted biomass pyramid, long seen as an exception, is increasingly recognized as widespread in certain systems—e.g., marine ecosystems (Sandin and Zgliczynski 2015; McCauley et al. 2018; Woodson, Schramski, and Joye 2018).
To investigate biomass distribution in our model, we use the Consumer to Resource biomass ratio, calculated as:
\(\displaystyle C{:}R = \frac{{C^{\ast}_i}}{{P^{\ast}_i}}\)
where the asterisk \(^{\ast}\) indicates the equilibrium biomass estimates for Consumer and Primary Producers and \(i \in [1, 2]\) indicates the ecosystem— either donor (1) or recipient (2). Values of \(C{:}R < 1\) identify ecosystems with a bottom-heavy biomass pyramid, whereas when \(C{:}R > 1\) the biomass pyramid is inverted and top-heavy (McCauley et al. 2018).
We begin by calculating the \(C{:}R\) ratio for each scenario investigated above.
First, we create a new dataframe bmDistr, where we will store the required
data: fertility conditions, movement rates values, efficiency rates values,
state variable values at equilibrium, and ecosystem function values. We will
also filter out those parameter sets that have no biological sense and are
unstable as per the stability analyses above.
Let’s look at the density distribution of these data. First, we exclude parameter sets that produce extreme values of \(C{:}R\)—that is, \(C{:}R > 10\).
Now, we pivot the resulting dataset and plot the density distribution of \(C{:}R\) values.
fert.labs <- c("Against-gradient", "Along-gradient", "Gradient-neutral")
names(fert.labs) <- c("Against", "Along", "Equal")
scale.labs <- c("Ecosystem 1", "Ecosystem 2", "Meta-ecosystem")
names(scale.labs) <- c("Ecosystem 1", "Ecosystem 2", "Meta-ecosystem")
ecoCR_long <- pivot_longer(ecoCR, cols = c(Eco1, Eco2, MetaEco), values_to = "CRratio",
names_to = "Scale") %>%
mutate_at(., vars(Scale), factor) %>%
mutate(., Scale = ifelse(Scale == "Eco1", "Ecosystem 1", ifelse(Scale ==
"Eco2", "Ecosystem 2", "Meta-ecosystem")))
ggplot(ecoCR_long, aes(x = CRratio, fill = Fertility, col = Fertility)) +
geom_density(aes(y = after_stat(density)), alpha = 0.1) + scale_color_highcontrast(reverse = T,
name = "Consumer movement", labels = c("Against-gradient", "Along-gradient",
"Gradient-neutral")) + scale_fill_highcontrast(reverse = T, name = "Consumer movement",
labels = c("Against-gradient", "Along-gradient", "Gradient-neutral")) +
facet_grid(Fertility ~ Scale, scales = "free", labeller = labeller(Fertility = fert.labs,
Scale = scale.labs)) + xlab("C:R ratio") + ylab("Density") + theme(legend.position = "top")
Figure 12: Density distribution when dataset is limited to C:R < 10. Note the different scales on the y-axis.
# ggsave(filename = '../Results/Density_CRratio_lessThan10.pdf',
# device = 'pdf', dpi = 300, width = 10, height = 7)
Here we look at how the \(C{:}R\) values change with increasing movement from
ecosystem 1, the donor, towards the dispersers’ pool Q (i.e., parameter g in
the model). The graph below uses the full dataset created above, ecoCR_long,
and fits linear models to each combination of Fertility scenario
(gradient-neutral, along-, and against-gradient) and Scale (ecosystem 1,
ecosystem 2, and meta-ecosystem).
ecoCRg_lm <- ecoCR_long %>%
nest_by(Fertility, Scale) %>%
mutate(model = list(lm(CRratio ~ g, data = data)))
ecoCRg_lmOutput <- ecoCRg_lm %>%
summarise(glance(model))
ecoCRg_lmOutput$g_slope <- NA
ecoCRg_lmOutput$g_se <- NA
for (i in 1:length(ecoCRg_lm[["model"]])) {
ecoCRg_lmOutput$g_slope[i] <- ecoCRg_lm$model[[i]]$coefficient[2]
ecoCRg_lmOutput$g_se[i] <- summary(ecoCRg_lm$model[[i]])$coefficient["g",
"Std. Error"]
}
ecoCRg_lmOutput %>%
gt(groupname_col = "Scale") %>%
fmt_number(columns = c(r.squared:p.value, logLik:deviance, g_slope:g_se),
decimals = 3) %>%
cols_label(., r.squared = "R^2", adj.r.squared = "Adj. R^2", sigma = "SE",
df.residual = "Residual df", nobs = "n", g_slope = "g Coeff.",
g_se = "g SE") %>%
gtsave(., filename = "../Results/CR_g_lm_res.rtf")
ecoCR_long %>%
ggplot(., aes(x = g, y = CRratio)) + geom_bin2d(bins = 70) + scale_fill_met_c("Hokusai2",
direction = 1) + geom_smooth(method = "lm", se = FALSE, col = "black",
linetype = 2, show.legend = T, size = 0.5) + geom_quantile(method = "rq",
quantiles = 0.5, col = "black", show.legend = T) + facet_grid(Fertility ~
Scale, scales = "free", labeller = labeller(Fertility = fert.labs,
Scale = scale.labs)) + ylab("C:R ratio") + xlab("Movement rate from Ecosystem 1 to Dispersers' Pool (g)")
Figure 13: Distribution of C:R values lower than 10 as the movement rate from ecosystem 1 to the dispersers’ pool (g) increases. Darker shades of blue correspond to a higher count of C:R value in a given bin. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data. C:R values appear to follow a decreasing trend as g increases in Ecosystem 1 and in the Meta-ecosystem, whereas they appear stable or slightly increasing in Ecosystem 2.
# ggtitle('Binned distribution of C:R', subtitle = 'C:R values
# [0,10].')
# ggsave(filename = '../Results/CRratioVg_binned.pdf', device =
# 'pdf', dpi = 300, width = 12, height = 12)
Now, let’s explore the relationships in a bit more detail. We do this by splitting the dataframe in two portions: \(C{:}R > 1\) and \(C{:}R < 1\). Here, we take a look at the relationship for \(C{:}R\) values \(> 1\) but \(< 10\).
ecoCRg_lm_above1 <- ecoCR_long %>%
filter(., CRratio > 1) %>%
nest_by(Fertility, Scale) %>%
mutate(model = list(lm(CRratio ~ g, data = data)))
ecoCRg_lmOutput_above1 <- ecoCRg_lm_above1 %>%
summarise(glance(model))
ecoCRg_lmOutput_above1$g_slope <- NA
ecoCRg_lmOutput_above1$g_se <- NA
for (i in 1:length(ecoCRg_lm_above1[["model"]])) {
ecoCRg_lmOutput_above1$g_slope[i] <- ecoCRg_lm_above1$model[[i]]$coefficient[2]
ecoCRg_lmOutput_above1$g_se[i] <- summary(ecoCRg_lm_above1$model[[i]])$coefficient["g",
"Std. Error"]
}
ecoCRg_lmOutput_above1 %>%
gt(groupname_col = "Scale") %>%
fmt_number(columns = c(r.squared:p.value, logLik:deviance, g_slope:g_se),
decimals = 3) %>%
cols_label(., r.squared = "R^2", adj.r.squared = "Adj. R^2", sigma = "SE",
df.residual = "Residual df", nobs = "n", g_slope = "g Coeff.",
g_se = "g SE") %>%
tab_header(title = md("Modeling C:R ratio vs Movement rate to Dispersers'
Pool (g)"),
subtitle = md("C:R values between (1,10].")) %>%
gtsave(., filename = "../Results/CR_g_lm_above1.rtf")
ecoCR_long %>%
filter(., CRratio > 1) %>%
ggplot(., aes(x = m, y = CRratio)) + geom_bin2d(bins = 70) + scale_fill_met_c("Hokusai2",
direction = 1) + geom_smooth(method = "lm", se = FALSE, col = "black",
linetype = 2, show.legend = T, size = 0.5) + geom_quantile(method = "rq",
quantiles = 0.5, col = "black", show.legend = T) + facet_grid(Fertility ~
Scale, scales = "free", labeller = labeller(Fertility = fert.labs,
Scale = scale.labs)) + ylab("C:R ratio") + xlab("Movement rate from Ecosystem 1 to Dispersers' Pool (g)")
Figure 14: Distribution of C:R values comprised between (1, 10] as the movement rate from ecosystem 1 to the dispersers’ pool (g) increases. As g increases, we do not observe any trend in the values of C:R. Note that C:R > 1 are generally indicative of top-heavy, inverted biomass pyramids. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.
# ggsave(filename = '../Results/CRratioVg_binned_above1.pdf', device
# = 'pdf', dpi = 300, width = 12, height = 12)
And here is the relationship for \(C{:}R\) values \(< 1\) but \(> 0\).
ecoCRg_lm_below1 <- ecoCR_long %>%
filter(., CRratio <= 1) %>%
nest_by(Fertility, Scale) %>%
mutate(model = list(lm(CRratio ~ g, data = data)))
ecoCRg_lmOutput_below1 <- ecoCRg_lm_below1 %>% summarise(glance(model))
ecoCRg_lmOutput_below1$g_slope <- NA
ecoCRg_lmOutput_below1$g_se <- NA
for (i in 1:length(ecoCRg_lm_below1[["model"]])) {
ecoCRg_lmOutput_below1$g_slope[i] <- ecoCRg_lm_below1$model[[i]]$coefficient[2]
ecoCRg_lmOutput_below1$g_se[i] <- summary(ecoCRg_lm_above1$model[[i]])$coefficient["g", "Std. Error"]
}
ecoCRg_lmOutput_below1 %>%
gt(groupname_col = "Scale") %>%
fmt_number(., columns = c(r.squared:p.value,logLik:deviance, g_slope:g_se),
decimals = 3) %>%
cols_label(., r.squared = "R^2",
adj.r.squared = "Adj. R^2",
sigma = "SE",
df.residual = "Residual df",
nobs = "n",
g_slope = "g Coeff.",
g_se = "g SE") %>%
tab_header(title = md("Modeling C:R ratio vs Movement rate to Dispersers'
Pool (g)"),
subtitle = md("C:R values between [0,1].")) %>%
gtsave(., filename = "../Results/CR_g_lm_below1.rtf")
ecoCR_long %>%
filter(., CRratio <= 1) %>%
ggplot(., aes(x = m, y = CRratio)) +
geom_bin2d(bins = 70) +
scale_fill_met_c("Hokusai2", direction = 1) +
geom_smooth(method = "lm", se = FALSE, col = "black",
linetype = 2, show.legend = T, size = .5) +
geom_quantile(method = "rq", quantiles = .5, col = "black",
show.legend = T) +
facet_grid(Fertility~Scale, scales = "free",
labeller = labeller(Fertility = fert.labs,
Scale = scale.labs)) +
ylab("C:R ratio") +
xlab("Movement rate from Ecosystem 1 to Dispersers' Pool (g)")
Figure 15: Distribution of C:R values [0,1] when movement rate from ecosystem 1 to the dispersers’ pool (g) increases. Again, no trend is apparent save for a seeming increase in C:R values in ecosystem 2 when consumers move along-gradient (central panel). Note, however, that this trend does not scale up to the meta-ecosystem. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.
#ggsave(filename = "../Results/CRratioVg_binned_below1.pdf", device = "pdf", dpi = 300, width = 12, height = 12)
Here, we investigate how the \(C{:}R\) biomass ratio changes
with increasing values of m, the movement rate of consumers leaving the
dispersers’ pool Q to enter ecosystem 2. Here is the general relationship.
ecoCRm_lm <- ecoCR_long %>%
nest_by(Fertility, Scale) %>%
mutate(model = list(lm(CRratio ~ m, data = data)))
ecoCRm_lmOutput <- ecoCRm_lm %>%
summarise(glance(model))
ecoCRm_lmOutput$m_slope <- NA
ecoCRm_lmOutput$m_se <- NA
for (i in 1:length(ecoCRm_lm[["model"]])) {
ecoCRm_lmOutput$m_slope[i] <- ecoCRm_lm$model[[i]]$coefficient[2]
ecoCRm_lmOutput$m_se[i] <- summary(ecoCRm_lm$model[[i]])$coefficient["m",
"Std. Error"]
}
ecoCRm_lmOutput %>%
gt(groupname_col = "Scale") %>%
fmt_number(columns = c(r.squared:p.value, logLik:deviance, m_slope:m_se),
decimals = 3) %>%
cols_label(., r.squared = "R^2", adj.r.squared = "Adj. R^2", sigma = "SE",
df.residual = "Residual df", nobs = "n", m_slope = "m Coeff.",
m_se = "m SE") %>%
gtsave(., filename = "../Results/CR_m_lm_res.rtf")
ecoCR_long %>%
ggplot(., aes(x = m, y = CRratio, )) + geom_bin2d(bins = 70) + scale_fill_met_c("Hokusai2",
direction = 1) + geom_smooth(method = "lm", se = FALSE, col = "black",
linetype = 2, show.legend = T, size = 0.5) + geom_quantile(method = "rq",
quantiles = 0.5, col = "black", show.legend = T) + facet_grid(Fertility ~
Scale, scales = "free", labeller = labeller(Fertility = fert.labs,
Scale = scale.labs)) + ylab("C:R ratio") + xlab("Movement rate from the Dispersers' Pool to Ecosystem 2 (m)")
Figure 16: Density distribution of C:R values [0,10] as consumer movement rate from the dispersers’ pool to ecosystem 2 (m) increases. No trend is clearly visible, aside from an apparent increase in C:R in ecosystem 2 when consumers move along-gradients (central panel). Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.
# ggsave(filename = '../Results/CRratioVm_binned.pdf', device =
# 'pdf', dpi = 300, width = 12, height = 12)
Here is the relationship for \(C{:}R\) values \(> 1\) but \(< 10\).
ecoCRm_lm_above1 <- ecoCR_long %>%
filter(., CRratio > 1) %>%
nest_by(Fertility, Scale) %>%
mutate(model = list(lm(CRratio ~ m, data = data)))
ecoCRm_lmOutput_above1 <- ecoCRm_lm_above1 %>%
summarise(glance(model))
ecoCRm_lmOutput_above1$m_slope <- NA
ecoCRm_lmOutput_above1$m_se <- NA
for (i in 1:length(ecoCRm_lm_above1[["model"]])) {
ecoCRm_lmOutput_above1$m_slope[i] <- ecoCRm_lm_above1$model[[i]]$coefficients[2]
ecoCRm_lmOutput_above1$m_se[i] <- summary(ecoCRm_lm_above1$model[[i]])$coefficient["m",
"Std. Error"]
}
ecoCRm_lmOutput_above1 %>%
gt(groupname_col = "Scale") %>%
fmt_number(columns = c(r.squared:p.value, logLik:deviance, m_slope:m_se),
decimals = 3) %>%
cols_label(., r.squared = "R^2", adj.r.squared = "Adj. R^2", sigma = "SE",
df.residual = "Residual df", nobs = "n", m_slope = "m Coeff.",
m_se = "m SE") %>%
tab_header(title = md("Modeling C:R ratio vs Movement rate from Dispersers'
Pool (m)"),
subtitle = md("C:R values between (1,10].")) %>%
gtsave(., filename = "../Results/CR_m_lm_above1.rtf")
ecoCR_long %>%
filter(., CRratio > 1) %>%
ggplot(., aes(x = m, y = CRratio)) + geom_bin2d(bins = 70) + scale_fill_met_c("Hokusai2",
direction = 1) + geom_smooth(method = "lm", se = FALSE, col = "black",
linetype = 2, show.legend = T, size = 0.5) + geom_quantile(method = "rq",
quantiles = 0.5, col = "black", show.legend = T) + facet_grid(Fertility ~
Scale, scales = "free", labeller = labeller(Fertility = fert.labs,
Scale = scale.labs)) + ylab("C:R ratio") + xlab("Movement rate from the Dispersers' Pool to Ecosystem 2 (m)")
Figure 17: Distribution of C:R (1, 10] as movement rate from the dispersers’ pool to ecosystem 2 (m) increases. Focusing on only C:R values indicative of inverted, top-heavy biomass pyramids shows weakly increasing trends in the value of C:R in ecosystem 1 when consumers move against-gradient (top-left panel), and in ecosystem 2 when consumers move along-gradient (central panel). We also note an apparent decrese in C:R values in ecosystem 1 as consumers move in an homogeneous meta-ecosystem (bottom-left panel). Note, however, that none of these trends appears to influence biomass dynamics at the meta-ecosystem scale (right-most column). Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.
# ggsave(filename = '../Results/CRratioVm_binned_above1.pdf', device
# = 'pdf', dpi = 300, width = 12, height = 12)
And here is the relationship for \(C{:}R\) values \(< 1\) but \(> 0\).
ecoCRm_lm_below1 <- ecoCR_long %>%
filter(., CRratio <= 1) %>%
nest_by(Fertility, Scale) %>%
mutate(model = list(lm(CRratio ~ m, data = data)))
ecoCRm_lmOutput_below1 <- ecoCRm_lm_below1 %>%
summarise(glance(model))
ecoCRm_lmOutput_below1$m_slope <- NA
ecoCRm_lmOutput_below1$m_se <- NA
for (i in 1:length(ecoCRm_lm_below1[["model"]])) {
ecoCRm_lmOutput_below1$m_slope[i] <- ecoCRm_lm_below1$model[[i]]$coefficient[2]
ecoCRm_lmOutput_below1$m_se[i] <- summary(ecoCRm_lm_below1$model[[i]])$coefficient["m",
"Std. Error"]
}
ecoCRm_lmOutput_below1 %>%
gt(groupname_col = "Scale") %>%
fmt_number(columns = c(r.squared:p.value, logLik:deviance, m_slope:m_se),
decimals = 3) %>%
cols_label(., r.squared = "R^2", adj.r.squared = "Adj. R^2", sigma = "SE",
df.residual = "Residual df", nobs = "n", m_slope = "m Coeff.",
m_se = "m SE") %>%
tab_header(title = md("Modeling C:R ratio vs Movement rate from Dispersers' Pool (m)"),
subtitle = md("C:R values between [0,1].")) %>%
gtsave(., filename = "../Results/CR_m_lm_below1.rtf")
ecoCR_long %>%
filter(., CRratio <= 1) %>%
ggplot(., aes(x = m, y = CRratio)) + geom_bin2d(bins = 70) + scale_fill_met_c("Hokusai2",
direction = 1) + geom_smooth(method = "lm", se = FALSE, col = "black",
linetype = 2, show.legend = T, size = 0.5) + geom_quantile(method = "rq",
quantiles = 0.5, col = "black", show.legend = T) + facet_grid(Fertility ~
Scale, scales = "free", labeller = labeller(Fertility = fert.labs,
Scale = scale.labs)) + ylab("C:R ratio") + xlab("Movement rate from the Dispersers' Pool to Ecosystem 2 (m)")
Figure 18: Distribution of C:R [0,1] as consumer movement rate from the dispersers pool to ecosystem 2 (m) increases. Here, we notice apparently increasing trends in ecosystem 2 when consumers move along-gradient (central panel) and in an homogeneous meta-ecosystem (center-bottom panel). We do not see any influence of these trends at the meta-ecosystem scale. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.
# ggsave(filename = '../Results/CRratioVm_binned_below1.pdf', device
# = 'pdf', dpi = 300, width = 12, height = 12)
Here, we investigate how the \(C{:}R\) biomass ratio changes
with increasing values of c, the death rate of consumers while in the
dispersers’ pool. Consumer biomass lost while transiting through Q is
lost from the meta-ecosystem and does not contribute to the recycling
pathways at either local or regional scale. Hence, it is important to
track how biomass distribution in the system may vary with this parameter. Here
is the general relationship.
ecoCRc_lm <- ecoCR_long %>%
nest_by(Fertility, Scale) %>%
mutate(model = list(lm(CRratio ~ c, data = data)))
ecoCRc_lmOutput <- ecoCRc_lm %>%
summarise(glance(model))
ecoCRc_lmOutput$c_slope <- NA
ecoCRc_lmOutput$c_se <- NA
for (i in 1:length(ecoCRc_lm[["model"]])) {
ecoCRc_lmOutput$c_slope[i] <- ecoCRc_lm$model[[i]]$coefficient[2]
ecoCRc_lmOutput$c_se[i] <- summary(ecoCRc_lm$model[[i]])$coefficient["c",
"Std. Error"]
}
ecoCRc_lmOutput %>%
gt(groupname_col = "Scale") %>%
fmt_number(columns = c(r.squared:p.value, logLik:deviance, c_slope:c_se),
decimals = 3) %>%
cols_label(., r.squared = "R^2", adj.r.squared = "Adj. R^2", sigma = "SE",
df.residual = "Residual df", nobs = "n", c_slope = "c Coeff.",
c_se = "c SE") %>%
gtsave(., filename = "../Results/CR_c_lm_res.rtf")
ecoCR_long %>%
ggplot(., aes(x = c, y = CRratio, )) + geom_bin2d(bins = 70) + scale_fill_met_c("Hokusai2",
direction = 1) + geom_smooth(method = "lm", se = FALSE, col = "black",
linetype = 2, show.legend = T, size = 0.5) + geom_quantile(method = "rq",
quantiles = 0.5, col = "black", show.legend = T) + facet_grid(Fertility ~
Scale, scales = "free", labeller = labeller(Fertility = fert.labs,
Scale = scale.labs)) + ylab("C:R ratio") + xlab("Consumers death rate while in the Dispersers' Pool (c)")
Figure 19: Distribution of C:R values [0,10] as the consumer death rate in the dispersers’ pool (c) increases. No trend is apparent, aside from a weak decrease in ecosystem 2 (central and center-bottom panels) that have no influence on meta-ecosystem biomass dynamics. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.
# ggsave(filename = '../Results/CRratioVc_binned.pdf', device =
# 'pdf', dpi = 300, width = 12, height = 12)
Here is the relationship for \(C{:}R\) values \(> 1\) but \(< 10\).
ecoCRc_lm_above1 <- ecoCR_long %>%
filter(., CRratio > 1) %>%
nest_by(Fertility, Scale) %>%
mutate(model = list(lm(CRratio ~ c, data = data)))
ecoCRc_lmOutput_above1 <- ecoCRc_lm_above1 %>%
summarise(glance(model))
ecoCRc_lmOutput_above1$c_slope <- NA
ecoCRc_lmOutput_above1$c_se <- NA
for (i in 1:length(ecoCRc_lm_above1[["model"]])) {
ecoCRc_lmOutput_above1$c_slope[i] <- ecoCRc_lm_above1$model[[i]]$coefficients[2]
ecoCRc_lmOutput_above1$c_se[i] <- summary(ecoCRc_lm_above1$model[[i]])$coefficient["c",
"Std. Error"]
}
ecoCRc_lmOutput_above1 %>%
gt(groupname_col = "Scale") %>%
fmt_number(columns = c(r.squared:p.value, logLik:deviance, c_slope:c_se),
decimals = 3) %>%
cols_label(., r.squared = "R^2", adj.r.squared = "Adj. R^2", sigma = "SE",
df.residual = "Residual df", nobs = "n", c_slope = "c Coeff.",
c_se = "c SE") %>%
tab_header(title = md("Modeling C:R ratio vs death rate while in the
Dispersers' Pool (c)"),
subtitle = md("C:R values between (1,10].")) %>%
gtsave(., filename = "../Results/CR_c_lm_above1.rtf")
ecoCR_long %>%
filter(., CRratio > 1) %>%
ggplot(., aes(x = c, y = CRratio)) + geom_bin2d(bins = 70) + scale_fill_met_c("Hokusai2",
direction = 1) + geom_smooth(method = "lm", se = FALSE, col = "black",
linetype = 2, show.legend = T, size = 0.5) + geom_quantile(method = "rq",
quantiles = 0.5, col = "black", show.legend = T) + facet_grid(Fertility ~
Scale, scales = "free", labeller = labeller(Fertility = fert.labs,
Scale = scale.labs)) + ylab("C:R ratio") + xlab("Consumer death rate while in the Dispersers' Pool (c)")
Figure 20: Distribution of C:R (1,10] as consumer death rate in the dispersers’ pool (c) increases. Regardless of the type of consumer movement, we see no trend in the distribution of C:R values that identify top-heavy, inverted biomass pyramids. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.
# ggsave(filename = '../Results/CRratioVc_binned_above1.pdf', device
# = 'pdf', dpi = 300, width = 12, height = 12)
And here is the relationship for \(C{:}R\) values \(< 1\) but \(> 0\).
ecoCRc_lm_below1 <- ecoCR_long %>%
filter(., CRratio <= 1) %>%
nest_by(Fertility, Scale) %>%
mutate(model = list(lm(CRratio ~ c, data = data)))
ecoCRc_lmOutput_below1 <- ecoCRc_lm_below1 %>%
summarise(glance(model))
ecoCRc_lmOutput_below1$c_slope <- NA
ecoCRc_lmOutput_below1$c_se <- NA
for (i in 1:length(ecoCRc_lm_below1[["model"]])) {
ecoCRc_lmOutput_below1$c_slope[i] <- ecoCRc_lm_below1$model[[i]]$coefficient[2]
ecoCRc_lmOutput_below1$c_se[i] <- summary(ecoCRc_lm_below1$model[[i]])$coefficient["c",
"Std. Error"]
}
ecoCRc_lmOutput_below1 %>%
gt(groupname_col = "Scale") %>%
fmt_number(columns = c(r.squared:p.value, logLik:deviance, c_slope:c_se),
decimals = 3) %>%
cols_label(., r.squared = "R^2", adj.r.squared = "Adj. R^2", sigma = "SE",
df.residual = "Residual df", nobs = "n", c_slope = "c Coeff.",
c_se = "c SE") %>%
tab_header(title = md("Modeling C:R ratio vs death rate while in the
Dispersers' pool (c)"),
subtitle = md("C:R values between [0,1].")) %>%
gtsave(., filename = "../Results/CR_c_lm_below1.rtf")
ecoCR_long %>%
filter(., CRratio <= 1) %>%
ggplot(., aes(x = c, y = CRratio)) + geom_bin2d(bins = 70) + scale_fill_met_c("Hokusai2",
direction = 1) + geom_smooth(method = "lm", se = FALSE, col = "black",
linetype = 2, show.legend = T, size = 0.5) + geom_quantile(method = "rq",
quantiles = 0.5, col = "black", show.legend = T) + facet_grid(Fertility ~
Scale, scales = "free", labeller = labeller(Fertility = fert.labs,
Scale = scale.labs)) + ylab("C:R ratio") + xlab("Consumer death rate while in the Dispersers' Pool (c)")
Figure 21: Distribution of C:R [0,1] as consumer death rate in the dispersers’ pool (c) increases. Focusing on C:R < 1 brings out a weak, negative trend in ecosystem 2 when consumers move along-gradient (central panel) and in an homogeneous meta-ecosystem (center-bottom panel). Neither of these trends have any influence at the meta-ecosystem scale. Lines correspond to a linear regression (black dashed line) and to a quantile regression calculated for the median (black continuous line) fitted to the data.
# ggsave(filename = '../Results/CRratioVc_binned_below1.pdf', device
# = 'pdf', dpi = 300, width = 12, height = 12)
Here we produce graphs of the biomass pyramids produced by consumer movement across ecosystems. These graphs allow for easy visualization of the changes in the distribution of biomass in either local ecosystem, or in the meta-ecosystem as a whole, following different types of consumer movement.
We produce four graphs, one each for the following four cases: (i) classic, bottom-heavy pyramid in the donor ecosystem 1; (ii) classic, bottom-heavy pyramid in the recipient ecosystem 2; (iii) inverted, top-heavy pyramid in the donor ecosystem 1; (iv) inverted, top-heavy pyramid in the recipient ecosystem 2. We separate between bottom- and top-heavy pyramids in either the donor or recipient ecosystems as these local conditions can influence the effects of consumer movement on biomass distribution and thus ecosystem functioning. Note that these four biomass-informed categories will return when we perform a global sensitivity analyses on our results below (Sensitivity Analyses).
The code we use here was initially developed by Dr. A. Z. Andis Arietta.
# first, produce a df that contains the values that will be printed
# on top of each bar in the biomass pyramids
summary_df <- ecoCR %>%
ungroup() %>%
select(biosense:I2, N1:Ctot, Eco1) %>%
filter(., Eco1 < 1) %>%
pivot_longer(cols = c(N1:Ctot), names_to = "Compartment", values_to = "Biomass") %>%
mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
"P1" | Compartment == "C1", "Ecosystem 1", ifelse(Compartment ==
"N2" | Compartment == "P2" | Compartment == "C2", "Ecosystem 2",
"Meta-ecosystem")), .before = Compartment, Compartment = ifelse(Compartment ==
"N1" | Compartment == "N2" | Compartment == "Ntot", "Nutrients",
ifelse(Compartment == "P1" | Compartment == "P2" | Compartment ==
"Ptot", "Primary\nProducers", "Consumers")), Compartment = fct_relevel(Compartment,
"Nutrients", "Primary\nProducers", "Consumers")) %>%
group_by(Fertility, Ecosystem, Compartment)
# then, produce the plot
ecoCR %>%
ungroup() %>%
select(biosense:I2, N1:Ctot, Eco1) %>%
filter(Eco1 < 1) %>%
pivot_longer(cols = c(N1:Ctot), names_to = "Compartment", values_to = "Biomass") %>%
mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
"P1" | Compartment == "C1", "Ecosystem 1", ifelse(Compartment ==
"N2" | Compartment == "P2" | Compartment == "C2", "Ecosystem 2",
"Meta-ecosystem")), .before = Compartment, Compartment = ifelse(Compartment ==
"N1" | Compartment == "N2" | Compartment == "Ntot", "Nutrients",
ifelse(Compartment == "P1" | Compartment == "P2" | Compartment ==
"Ptot", "Primary\nProducers", "Consumers")), Compartment = fct_relevel(Compartment,
"Nutrients", "Primary\nProducers", "Consumers")) %>%
ggplot(aes(x = Compartment, y = Biomass)) + geom_bar(aes(fill = Compartment,
color = Compartment), stat = "identity", width = 0.9) + geom_bar(data = . %>%
mutate(Biomass = -Biomass), aes(fill = Compartment, color = Compartment),
stat = "identity", width = 0.9) + stat_summary(data = summary_df, geom = "label",
fun.args = list(mult = 1, na.rm = T), aes(label = round(after_stat(y),
2)), position = position_dodge(width = 1), show.legend = F, size = 2,
label.padding = unit(0.1, "lines"), label.size = 0.1, color = "black") +
scale_color_met_d(name = "Isfahan2", direction = 1) + scale_fill_met_d(name = "Isfahan2",
direction = 1) + facet_grid(Fertility ~ Ecosystem, scales = "free") +
coord_flip() + theme_classic() + theme(text = element_text(size = 14),
axis.title = element_blank(), axis.ticks = element_blank(), axis.text = element_blank(),
legend.position = "top")
Figure 22: Biomass pyramids for ecosystem 1 (donor), ecosystem 2 (recipient), and the meta-ecosystem when the Consumer:Resource ratio in ecosystem 1 is < 1. Bottom-heavy biomass pyramids dominate the system in all movement scenarios. However, note how against-gradient consumer movement consistently allows for—on average—higher levels of Primary Producers biomass compared to both along-gradient and gradient-neutral consumer movement. The labels inside the bars show the average biomass for each ecosystem compartment.
ggsave(filename = "../Results/Bpyram_CRbelow1_eco1.pdf", dpi = 300, device = "pdf")
# first, produce a df that contains the values that will be printed
# on top of each bar in the biomass pyramids
summary_df <- ecoCR %>%
ungroup() %>%
select(biosense:I2, N1:Ctot, Eco2) %>%
filter(Eco2 < 1) %>%
pivot_longer(cols = c(N1:Ctot), names_to = "Compartment", values_to = "Biomass") %>%
mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
"P1" | Compartment == "C1", "Ecosystem 1", ifelse(Compartment ==
"N2" | Compartment == "P2" | Compartment == "C2", "Ecosystem 2",
"Meta-ecosystem")), .before = Compartment, Compartment = ifelse(Compartment ==
"N1" | Compartment == "N2" | Compartment == "Ntot", "Nutrients",
ifelse(Compartment == "P1" | Compartment == "P2" | Compartment ==
"Ptot", "Primary\nProducers", "Consumers")), Compartment = fct_relevel(Compartment,
"Nutrients", "Primary\nProducers", "Consumers")) %>%
group_by(Fertility, Ecosystem, Compartment)
# then, produce the plot
ecoCR %>%
ungroup() %>%
select(biosense:I2, N1:Ctot, Eco2) %>%
filter(Eco2 < 1) %>%
pivot_longer(cols = c(N1:Ctot), names_to = "Compartment", values_to = "Biomass") %>%
mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
"P1" | Compartment == "C1", "Ecosystem 1", ifelse(Compartment ==
"N2" | Compartment == "P2" | Compartment == "C2", "Ecosystem 2",
"Meta-ecosystem")), .before = Compartment, Compartment = ifelse(Compartment ==
"N1" | Compartment == "N2" | Compartment == "Ntot", "Nutrients",
ifelse(Compartment == "P1" | Compartment == "P2" | Compartment ==
"Ptot", "Primary\nProducers", "Consumers")), Compartment = fct_relevel(Compartment,
"Nutrients", "Primary\nProducers", "Consumers")) %>%
ggplot(aes(x = Compartment, y = Biomass)) + geom_bar(aes(fill = Compartment,
color = Compartment), stat = "identity", width = 0.9) + geom_bar(data = . %>%
mutate(Biomass = -Biomass), aes(fill = Compartment, color = Compartment),
stat = "identity", width = 0.9) + stat_summary(data = summary_df, geom = "label",
fun.args = list(mult = 1, na.rm = T), aes(label = round(after_stat(y),
2)), position = position_dodge(width = 1), show.legend = F, size = 2,
label.padding = unit(0.1, "lines"), label.size = 0.1, color = "black") +
scale_color_met_d(name = "Isfahan2", direction = 1) + scale_fill_met_d(name = "Isfahan2",
direction = 1) + facet_grid(Fertility ~ Ecosystem, scales = "free") +
coord_flip() + theme_classic() + theme(text = element_text(size = 14),
axis.title = element_blank(), axis.ticks = element_blank(), axis.text = element_blank(),
legend.position = "top")
Figure 23: Biomass pyramids for ecosystem 1 (donor), ecosystem 2 (recipient), and the meta-ecosystem when the Consumer:Resource ratio in ecosystem 2 is < 1. In this case, the biomass distribution in the system is more nuanced and variable, and depends on the direction of consumer movement. Against-gradient consumer movement leads to comparable average biomass of Primary Producers in the two ecosystems. Conversely, in the along-gradient and gradient-neutral scenarios, the distribution of Primary Producers’ biomass is skewed towards ecosystem 1 (the donor). Note also how the distribution of Nutrients varies among donor and recipient ecosystems across movement scenarios.
ggsave(filename = "../Results/Bpyram_CRbelow1_eco2.pdf", dpi = 300, device = "pdf")
# first, produce a df that contains the values that will be printed
# on top of each bar in the biomass pyramids
summary_df <- ecoCR %>%
ungroup() %>%
select(biosense:I2, N1:Ctot, Eco1) %>%
filter(Eco1 > 1) %>%
pivot_longer(cols = c(N1:Ctot), names_to = "Compartment", values_to = "Biomass") %>%
mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
"P1" | Compartment == "C1", "Ecosystem 1", ifelse(Compartment ==
"N2" | Compartment == "P2" | Compartment == "C2", "Ecosystem 2",
"Meta-ecosystem")), .before = Compartment, Compartment = ifelse(Compartment ==
"N1" | Compartment == "N2" | Compartment == "Ntot", "Nutrients",
ifelse(Compartment == "P1" | Compartment == "P2" | Compartment ==
"Ptot", "Primary\nProducers", "Consumers")), Compartment = fct_relevel(Compartment,
"Nutrients", "Primary\nProducers", "Consumers")) %>%
group_by(Fertility, Ecosystem, Compartment)
# then, produce the plot
ecoCR %>%
ungroup() %>%
select(biosense:I2, N1:Ctot, Eco1) %>%
filter(Eco1 > 1) %>%
pivot_longer(cols = c(N1:Ctot), names_to = "Compartment", values_to = "Biomass") %>%
mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
"P1" | Compartment == "C1", "Ecosystem 1", ifelse(Compartment ==
"N2" | Compartment == "P2" | Compartment == "C2", "Ecosystem 2",
"Meta-ecosystem")), .before = Compartment, Compartment = ifelse(Compartment ==
"N1" | Compartment == "N2" | Compartment == "Ntot", "Nutrients",
ifelse(Compartment == "P1" | Compartment == "P2" | Compartment ==
"Ptot", "Primary\nProducers", "Consumers")), Compartment = fct_relevel(Compartment,
"Nutrients", "Primary\nProducers", "Consumers")) %>%
ggplot(aes(x = Compartment, y = Biomass)) + geom_bar(aes(fill = Compartment,
color = Compartment), stat = "identity", width = 0.9) + geom_bar(data = . %>%
mutate(Biomass = -Biomass), aes(fill = Compartment, color = Compartment),
stat = "identity", width = 0.9) + stat_summary(data = summary_df, geom = "label",
fun.args = list(mult = 1, na.rm = T), aes(label = round(after_stat(y),
2)), position = position_dodge(width = 1), show.legend = F, size = 2,
label.padding = unit(0.1, "lines"), label.size = 0.1, color = "black") +
scale_color_met_d(name = "Isfahan2", direction = 1) + scale_fill_met_d(name = "Isfahan2",
direction = 1) + facet_grid(Fertility ~ Ecosystem, scales = "free") +
coord_flip() + theme_classic() + theme(text = element_text(size = 14),
axis.title = element_blank(), axis.ticks = element_blank(), axis.text = element_blank(),
legend.position = "top")
Figure 24: Biomass pyramids for ecosystem 1 (donor), ecosystem 2 (recipient), and the meta-ecosystem when the Consumer:Resource ratio in ecosystem 1 is > 1. In this case, while we see top-heavy, inverted biomass pyramids in the donor ecosystem, bottom-heavy pyramids still dominate the recipient ecosystem and ultimately the meta-ecosystem as a whole. In the against-gradient scenario, we can see an accumulation of nutrients in the recipient ecosystem (resource pooling effect), whereas the along-gradient scenario leads to a more equal distribution of nutrients in the system that resembles what we find in the gradient-neutral scenario.
ggsave(filename = "../Results/Bpyram_CRabove1_eco1.pdf", dpi = 300, device = "pdf")
# first, produce a df that contains the values that will be printed
# on top of each bar in the biomass pyramids
summary_df <- ecoCR %>%
ungroup() %>%
select(biosense:I2, N1:Ctot, Eco2) %>%
filter(Eco2 > 1) %>%
pivot_longer(cols = c(N1:Ctot), names_to = "Compartment", values_to = "Biomass") %>%
mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
"P1" | Compartment == "C1", "Ecosystem 1", ifelse(Compartment ==
"N2" | Compartment == "P2" | Compartment == "C2", "Ecosystem 2",
"Meta-ecosystem")), .before = Compartment, Compartment = ifelse(Compartment ==
"N1" | Compartment == "N2" | Compartment == "Ntot", "Nutrients",
ifelse(Compartment == "P1" | Compartment == "P2" | Compartment ==
"Ptot", "Primary\nProducers", "Consumers")), Compartment = fct_relevel(Compartment,
"Nutrients", "Primary\nProducers", "Consumers")) %>%
group_by(Fertility, Ecosystem, Compartment)
# then, produce the plot
ecoCR %>%
ungroup() %>%
select(biosense:I2, N1:Ctot, Eco2) %>%
filter(Eco2 > 1) %>%
pivot_longer(cols = c(N1:Ctot), names_to = "Compartment", values_to = "Biomass") %>%
mutate(., Ecosystem = if_else(Compartment == "N1" | Compartment ==
"P1" | Compartment == "C1", "Ecosystem 1", ifelse(Compartment ==
"N2" | Compartment == "P2" | Compartment == "C2", "Ecosystem 2",
"Meta-ecosystem")), .before = Compartment, Compartment = ifelse(Compartment ==
"N1" | Compartment == "N2" | Compartment == "Ntot", "Nutrients",
ifelse(Compartment == "P1" | Compartment == "P2" | Compartment ==
"Ptot", "Primary\nProducers", "Consumers")), Compartment = fct_relevel(Compartment,
"Nutrients", "Primary\nProducers", "Consumers")) %>%
ggplot(aes(x = Compartment, y = Biomass)) + geom_bar(aes(fill = Compartment,
color = Compartment), stat = "identity", width = 0.9) + geom_bar(data = . %>%
mutate(Biomass = -Biomass), aes(fill = Compartment, color = Compartment),
stat = "identity", width = 0.9) + stat_summary(data = summary_df, geom = "label",
fun.args = list(mult = 1, na.rm = T), aes(label = round(after_stat(y),
2)), position = position_dodge(width = 1), show.legend = F, size = 2,
label.padding = unit(0.1, "lines"), label.size = 0.1, color = "black") +
scale_color_met_d(name = "Isfahan2", direction = 1) + scale_fill_met_d(name = "Isfahan2",
direction = 1) + facet_grid(Fertility ~ Ecosystem, scales = "free") +
coord_flip() + theme_classic() + theme(text = element_text(size = 14),
axis.title = element_blank(), axis.ticks = element_blank(), axis.text = element_blank(),
legend.position = "top")
Figure 25: Biomass pyramids for ecosystem 1 (donor), ecosystem 2 (recipient), and the meta-ecosystem when the Consumer:Resource ratio in ecosystem 2 is > 1. Here, we observe the spatial trophic cascade in action as a lack of top-down control on Primary Producers leads to a highly skewed accumulation of biomass for this compartment in ecosystem 1, the donor, whereas ecosystem 2 shows much reduced autothroph presence. Conversely, both Consumer biomass and Nutrient stock are higher in ecosystem 2, the recipient, in all movement scenarios. Notice how, when the recipient ecosystem shows a top-heavy biomass pyramid, along-gradient consumer movement is the only scenario in which the distribution of nutrients is somewhat homogeneous across the landscape.
ggsave(filename = "../Results/Bpyram_CRabove1_eco2.pdf", dpi = 300, device = "pdf")
The three tabs below contain tables showing summary statistics for the linear models fitted to the distribution of biomass in the meta-ecosystem. The table shown in tab [Movement from ecosystem 1] is intended as a companion to Figures 13, 14, and 15. Similarly, the table in tab [Movement towards ecosystem 2] is a companion to Figures 16, 17, and 18. Finally, the table in tab Death rate in Dispersers’ pool provides details to accompany Figures 19, 20, and 21.
# first, compute median of log10 ratio values
ecoCR_g_table <- bind_rows(`above 1` = ecoCRg_lmOutput_above1, `below 1` = ecoCRg_lmOutput_below1,
.id = "CRslice") %>%
ungroup() %>%
select(., CRslice:adj.r.squared, g_slope:g_se) %>%
pivot_wider(., names_from = "CRslice", values_from = c(r.squared, adj.r.squared,
g_slope, g_se)) %>%
gt(., groupname_col = "Scale") %>%
tab_spanner(label = md("**C:R (1, 10]**"), columns = c("r.squared_above 1",
"adj.r.squared_above 1", "g_slope_above 1", "g_se_above 1")) %>%
tab_spanner(label = md("**C:R [0, 1]**"), columns = c("r.squared_below 1",
"adj.r.squared_below 1", "g_slope_below 1", "g_se_below 1")) %>%
fmt_number(., columns = 3:10, decimals = 3) %>%
cols_label(., `r.squared_above 1` = md("R sq."), `adj.r.squared_above 1` = md("Adj. R sq."),
`g_slope_above 1` = md("g Coeff."), `g_se_above 1` = md("g SE"),
`r.squared_below 1` = md("R sq."), `adj.r.squared_below 1` = md("Adj. R sq."),
`g_slope_below 1` = md("g Coeff."), `g_se_below 1` = md("g SE")) %>%
# data_color(., columns = c(3:10), colors = col_numeric(palette =
# 'PiYG', domain = c(-1,1))) %>%
tab_header(title = md("**Relationships between biomass _C:R_ and movement from ecosystem 1.**"),
subtitle = md("The table shows the slope value (g Coeff.) and standard error (g SE) for linear regression models fitted to the biomass _C:R_ data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between _C:R_ values larger than 1 but below 10, and _C:R_ values below 1 but above 0.")) %>%
tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())
ecoCR_g_table
| Relationships between biomass C:R and movement from ecosystem 1. | ||||||||
| The table shows the slope value (g Coeff.) and standard error (g SE) for linear regression models fitted to the biomass C:R data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between C:R values larger than 1 but below 10, and C:R values below 1 but above 0. | ||||||||
| Fertility | C:R (1, 10] | C:R [0, 1] | ||||||
|---|---|---|---|---|---|---|---|---|
| R sq. | Adj. R sq. | g Coeff. | g SE | R sq. | Adj. R sq. | g Coeff. | g SE | |
| Ecosystem 1 | ||||||||
| Against | 0.083 | 0.078 | −0.728 | 0.180 | 0.119 | 0.119 | −0.017 | 0.180 |
| Along | 0.148 | 0.147 | −0.343 | 0.022 | 0.010 | 0.010 | −0.009 | 0.022 |
| Equal | 0.177 | 0.176 | −0.537 | 0.040 | 0.039 | 0.039 | −0.016 | 0.040 |
| Ecosystem 2 | ||||||||
| Against | 0.001 | 0.001 | −0.025 | 0.016 | 0.000 | 0.000 | 0.001 | 0.016 |
| Along | 0.000 | 0.000 | 0.018 | 0.015 | 0.006 | 0.005 | 0.007 | 0.015 |
| Equal | 0.001 | 0.001 | −0.027 | 0.015 | 0.001 | 0.000 | 0.003 | 0.015 |
| Meta-ecosystem | ||||||||
| Against | 0.034 | 0.033 | −0.078 | 0.012 | 0.002 | 0.001 | −0.004 | 0.012 |
| Along | 0.048 | 0.047 | −0.100 | 0.011 | 0.002 | 0.002 | −0.004 | 0.011 |
| Equal | 0.049 | 0.048 | −0.106 | 0.012 | 0.005 | 0.005 | −0.006 | 0.012 |
# first, compute median of log10 ratio values
ecoCR_m_table <- bind_rows(`above 1` = ecoCRm_lmOutput_above1, `below 1` = ecoCRm_lmOutput_below1,
.id = "CRslice") %>%
ungroup() %>%
select(., CRslice:adj.r.squared, m_slope:m_se) %>%
pivot_wider(., names_from = "CRslice", values_from = c(r.squared, adj.r.squared,
m_slope, m_se)) %>%
gt(., groupname_col = "Scale") %>%
tab_spanner(label = md("**C:R (1, 10]**"), columns = c("r.squared_above 1",
"adj.r.squared_above 1", "m_slope_above 1", "m_se_above 1")) %>%
tab_spanner(label = md("**C:R [0, 1]**"), columns = c("r.squared_below 1",
"adj.r.squared_below 1", "m_slope_below 1", "m_se_below 1")) %>%
fmt_number(., columns = 3:10, decimals = 3) %>%
cols_label(., `r.squared_above 1` = md("R sq."), `adj.r.squared_above 1` = md("Adj. R sq."),
`m_slope_above 1` = md("m Coeff."), `m_se_above 1` = md("m SE"),
`r.squared_below 1` = md("R sq."), `adj.r.squared_below 1` = md("Adj. R sq."),
`m_slope_below 1` = md("m Coeff."), `m_se_below 1` = md("m SE")) %>%
# data_color(., columns = c(3:10), colors = col_numeric(palette =
# 'PiYG', domain = c(-,))) %>%
tab_header(title = md("**Relationships between biomass _C:R_ and movement into ecosystem 2.**"),
subtitle = md("The table shows the slope value (m Coeff.) and standard error (m SE) for linear regression models fitted to the biomass _C:R_ data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between _C:R_ values larger than 1 but below 10, and _C:R_ values below 1 but above 0.")) %>%
tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())
ecoCR_m_table
| Relationships between biomass C:R and movement into ecosystem 2. | ||||||||
| The table shows the slope value (m Coeff.) and standard error (m SE) for linear regression models fitted to the biomass C:R data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between C:R values larger than 1 but below 10, and C:R values below 1 but above 0. | ||||||||
| Fertility | C:R (1, 10] | C:R [0, 1] | ||||||
|---|---|---|---|---|---|---|---|---|
| R sq. | Adj. R sq. | m Coeff. | m SE | R sq. | Adj. R sq. | m Coeff. | m SE | |
| Ecosystem 1 | ||||||||
| Against | 0.003 | −0.002 | 0.036 | 0.048 | 0.000 | 0.000 | 0.000 | 0.001 |
| Along | 0.000 | 0.000 | −0.014 | 0.018 | 0.001 | 0.001 | −0.002 | 0.001 |
| Equal | 0.002 | 0.001 | −0.031 | 0.026 | 0.000 | 0.000 | −0.001 | 0.001 |
| Ecosystem 2 | ||||||||
| Against | 0.000 | 0.000 | −0.008 | 0.016 | 0.000 | 0.000 | 0.000 | 0.002 |
| Along | 0.001 | 0.001 | 0.026 | 0.015 | 0.016 | 0.016 | 0.012 | 0.001 |
| Equal | 0.000 | 0.000 | −0.001 | 0.015 | 0.006 | 0.006 | 0.007 | 0.001 |
| Meta-ecosystem | ||||||||
| Against | 0.000 | −0.001 | 0.000 | 0.013 | 0.000 | 0.000 | 0.000 | 0.001 |
| Along | 0.000 | 0.000 | 0.010 | 0.011 | 0.000 | 0.000 | 0.001 | 0.001 |
| Equal | 0.001 | 0.000 | −0.014 | 0.013 | 0.000 | 0.000 | 0.001 | 0.001 |
# first, compute median of log10 ratio values
ecoCR_c_table <- bind_rows(`above 1` = ecoCRc_lmOutput_above1, `below 1` = ecoCRc_lmOutput_below1,
.id = "CRslice") %>%
ungroup() %>%
select(., CRslice:adj.r.squared, c_slope:c_se) %>%
pivot_wider(., names_from = "CRslice", values_from = c(r.squared, adj.r.squared,
c_slope, c_se)) %>%
gt(., groupname_col = "Scale") %>%
tab_spanner(label = md("**C:R (1, 10]**"), columns = c("r.squared_above 1",
"adj.r.squared_above 1", "c_slope_above 1", "c_se_above 1")) %>%
tab_spanner(label = md("**C:R [0, 1]**"), columns = c("r.squared_below 1",
"adj.r.squared_below 1", "c_slope_below 1", "c_se_below 1")) %>%
fmt_number(., columns = 3:10, decimals = 3) %>%
cols_label(., `r.squared_above 1` = md("R sq."), `adj.r.squared_above 1` = md("Adj. R sq."),
`c_slope_above 1` = md("c Coeff."), `c_se_above 1` = md("c SE"),
`r.squared_below 1` = md("R sq."), `adj.r.squared_below 1` = md("Adj. R sq."),
`c_slope_below 1` = md("c Coeff."), `c_se_below 1` = md("c SE")) %>%
# data_color(., columns = c(3:10), colors = col_numeric(palette =
# 'PiYG', domain = c(-,))) %>%
tab_header(title = md("**Relationships between biomass _C:R_ and death rate in dispersers' pool.**"),
subtitle = md("The table shows the slope value (c Coeff.) and standard error (c SE) for linear regression models fitted to the biomass _C:R_ data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between _C:R_ values larger than 1 but below 10, and _C:R_ values below 1 but above 0.")) %>%
tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())
ecoCR_c_table
| Relationships between biomass C:R and death rate in dispersers' pool. | ||||||||
| The table shows the slope value (c Coeff.) and standard error (c SE) for linear regression models fitted to the biomass C:R data. R squared and Adjusted R squared values provide a measure of fit of the model to the data. Data are split between C:R values larger than 1 but below 10, and C:R values below 1 but above 0. | ||||||||
| Fertility | C:R (1, 10] | C:R [0, 1] | ||||||
|---|---|---|---|---|---|---|---|---|
| R sq. | Adj. R sq. | c Coeff. | c SE | R sq. | Adj. R sq. | c Coeff. | c SE | |
| Ecosystem 1 | ||||||||
| Against | 0.005 | −0.001 | 0.047 | 0.050 | 0.000 | 0.000 | −0.001 | 0.001 |
| Along | 0.000 | 0.000 | −0.012 | 0.019 | 0.000 | 0.000 | 0.002 | 0.001 |
| Equal | 0.000 | −0.001 | 0.000 | 0.026 | 0.000 | 0.000 | 0.001 | 0.001 |
| Ecosystem 2 | ||||||||
| Against | 0.000 | 0.000 | 0.010 | 0.016 | 0.000 | 0.000 | 0.000 | 0.002 |
| Along | 0.000 | 0.000 | −0.004 | 0.014 | 0.002 | 0.002 | −0.004 | 0.001 |
| Equal | 0.001 | 0.001 | −0.030 | 0.015 | 0.002 | 0.002 | −0.005 | 0.002 |
| Meta-ecosystem | ||||||||
| Against | 0.003 | 0.003 | 0.025 | 0.012 | 0.000 | 0.000 | 0.000 | 0.001 |
| Along | 0.000 | −0.001 | −0.003 | 0.011 | 0.001 | 0.001 | −0.003 | 0.001 |
| Equal | 0.000 | −0.001 | 0.004 | 0.012 | 0.001 | 0.001 | −0.003 | 0.001 |
Here, we produce summary tables that show how the effects of consumer movement on local and meta-ecosystem functions vary as the biomass distribution changes in the system as a result of the consumer movement itself.
We begin by collating together the required dataframes.
Then, we lengthen the dataset as this is the preferred format for producing
tables with package gt.
# the first step in the pipe remove the by-ecosystem CR data, as
# these are now store in the CR column
CRcomparison_long <- select(CRcomparison, !c(biosense, stable, CR_Eco1,
CR_Eco2, CR_MetaEco)) %>%
pivot_longer(., cols = N1rr:PROD_Crr, names_to = "Compartment", values_to = "Ratio") %>%
dplyr::mutate(., Function = if_else(Compartment == "FLUX_P1rr" | Compartment ==
"FLUX_P2rr" | Compartment == "FLUX_C1rr" | Compartment == "FLUX_C2rr" |
Compartment == "FLUX_Prr" | Compartment == "FLUX_Crr", "Nutrient Flux",
if_else(Compartment == "PROD_P1rr" | Compartment == "PROD_C1rr" |
Compartment == "PROD_P2rr" | Compartment == "PROD_C2rr" | Compartment ==
"PROD_Prr" | Compartment == "PROD_Crr", "Trophic Productivity",
"Stock"))) %>%
dplyr::mutate(., Scale = if_else(Compartment == "N1rr" | Compartment ==
"P1rr" | Compartment == "C1rr" | Compartment == "FLUX_P1rr" | Compartment ==
"FLUX_C1rr" | Compartment == "PROD_P1rr" | Compartment == "PROD_C1rr",
"Ecosystem 1", if_else(Compartment == "N2rr" | Compartment == "P2rr" |
Compartment == "C2rr" | Compartment == "FLUX_P2rr" | Compartment ==
"FLUX_C2rr" | Compartment == "PROD_P2rr" | Compartment == "PROD_C2rr",
"Ecosystem 2", "Meta-ecosystem"))) %>%
dplyr::mutate(., Fertility = factor(Fertility), Function = factor(Function),
Scale = factor(Scale), Compartment = factor(Compartment))
We now remove all values of the Response Ratio that are \(< 0\) and \(log_{10}\)-transform them before calculating the median.
Finally, we produce summary tables that show the differences in median LRR when considering all stable parameter sets, parameter sets that produce \(C{:}R > 1\), and parameter sets that produce \(C{:}R < 1\). As before, we produce one table per ecosystem function.
CRsummStock <- CRcomparison_tables %>%
ungroup() %>%
filter(., Function == "Stock") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, Consumers = "C1rr",
Consumers = "C2rr", Consumers = "Crr", Producers = "P1rr", Producers = "P2rr",
Producers = "Prr", Nutrients = "N1rr", Nutrients = "N2rr", Nutrients = "Nrr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Nutrients",
"Producers", "Consumers")) %>%
pivot_wider(., names_from = c(CR, Fertility), values_from = median) %>%
gt(., groupname_col = "Scale") %>%
fmt_number(., columns = 4:9, decimals = 2) %>%
tab_spanner(label = md("**Along-gradient**"), columns = c("all_Along",
"above_1_Along", "below_1_Along")) %>%
tab_spanner(label = md("**Against-gradient**"), columns = c("all_Against",
"above_1_Against", "below_1_Against")) %>%
cols_label(Compartment = md("**Compartment**"), all_Against = md("*All C:R*"),
above_1_Against = md("*C:R > 1*"), below_1_Against = md("*C:R < 1*"),
all_Along = md("*All C:R*"), above_1_Along = md("*C:R > 1*"), below_1_Along = md("*C:R < 1*")) %>%
cols_hide("Function") %>%
tab_footnote(footnote = "Includes parameter sets producing C:R > 10.",
locations = cells_column_labels(columns = c("all_Against", "all_Along"))) %>%
tab_footnote(footnote = "Does not include parameter sets producing C:R > 10.",
locations = cells_column_labels(columns = c("above_1_Against",
"above_1_Along"))) %>%
tab_header(title = md("**Median LRR values of biomass and nutrient stocks**"),
subtitle = md("LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing _C:R_ > 1, or only parameter sets producing _C:R_ < 1.")) %>%
tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())
CRsummStock
| Median LRR values of biomass and nutrient stocks | ||||||
| LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing C:R > 1, or only parameter sets producing C:R < 1. | ||||||
| Compartment | Against-gradient | Along-gradient | ||||
|---|---|---|---|---|---|---|
| All C:R1 | C:R > 12 | C:R < 1 | All C:R1 | C:R > 12 | C:R < 1 | |
| Ecosystem 1 | ||||||
| Consumers | −0.72 | −0.71 | −0.72 | 0.26 | 0.26 | 0.26 |
| Nutrients | −0.16 | −0.55 | −0.15 | 0.12 | 0.20 | 0.08 |
| Producers | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |
| Ecosystem 2 | ||||||
| Consumers | 0.17 | 0.15 | 0.10 | −0.27 | −0.13 | −0.36 |
| Nutrients | 0.11 | 0.12 | 0.02 | −0.14 | −0.09 | −0.10 |
| Producers | 0.04 | 0.03 | 0.10 | −0.12 | −0.12 | −0.15 |
| Meta-ecosystem | ||||||
| Consumers | 0.08 | −0.15 | −0.12 | −0.09 | 0.04 | −0.09 |
| Nutrients | 0.03 | −0.18 | −0.04 | −0.03 | 0.06 | 0.00 |
| Producers | 0.01 | 0.02 | 0.05 | −0.02 | −0.05 | −0.04 |
| 1 Includes parameter sets producing C:R > 10. | ||||||
| 2 Does not include parameter sets producing C:R > 10. | ||||||
CRsummFlux <- CRcomparison_tables %>%
ungroup() %>%
filter(., Function == "Nutrient Flux") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, Consumers = "FLUX_C1rr",
Consumers = "FLUX_C2rr", Consumers = "FLUX_Crr", Producers = "FLUX_P1rr",
Producers = "FLUX_P2rr", Producers = "FLUX_Prr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Producers",
"Consumers")) %>%
pivot_wider(., names_from = c(CR, Fertility), values_from = median) %>%
gt(., groupname_col = "Scale") %>%
fmt_number(., columns = 4:9, decimals = 2) %>%
tab_spanner(label = md("**Along-gradient**"), columns = c("all_Along",
"above_1_Along", "below_1_Along")) %>%
tab_spanner(label = md("**Against-gradient**"), columns = c("all_Against",
"above_1_Against", "below_1_Against")) %>%
cols_label(Compartment = md("**Compartment**"), all_Against = md("*All C:R*"),
above_1_Against = md("*C:R > 1*"), below_1_Against = md("*C:R < 1*"),
all_Along = md("*All C:R*"), above_1_Along = md("*C:R > 1*"), below_1_Along = md("*C:R < 1*")) %>%
cols_hide("Function") %>%
tab_footnote(footnote = "Includes parameter sets producing C:R > 10.",
locations = cells_column_labels(columns = c("all_Against", "all_Along"))) %>%
tab_footnote(footnote = "Does not include parameter sets producing C:R > 10.",
locations = cells_column_labels(columns = c("above_1_Against",
"above_1_Along"))) %>%
tab_header(title = md("**Median LRR values of nutrient flux**"), subtitle = md("LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing _C:R_ > 1, or only parameter sets producing _C:R_ < 1.")) %>%
tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())
CRsummFlux
| Median LRR values of nutrient flux | ||||||
| LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing C:R > 1, or only parameter sets producing C:R < 1. | ||||||
| Compartment | Against-gradient | Along-gradient | ||||
|---|---|---|---|---|---|---|
| All C:R1 | C:R > 12 | C:R < 1 | All C:R1 | C:R > 12 | C:R < 1 | |
| Ecosystem 1 | ||||||
| Consumers | −0.72 | −0.71 | −0.72 | 0.26 | 0.26 | 0.26 |
| Producers | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 |
| Ecosystem 2 | ||||||
| Consumers | 0.17 | 0.15 | 0.10 | −0.27 | −0.13 | −0.36 |
| Producers | 0.04 | 0.03 | 0.10 | −0.12 | −0.12 | −0.15 |
| Meta-ecosystem | ||||||
| Consumers | 0.07 | −0.01 | −0.07 | −0.09 | 0.02 | −0.11 |
| Producers | 0.01 | 0.02 | 0.04 | −0.02 | −0.05 | −0.04 |
| 1 Includes parameter sets producing C:R > 10. | ||||||
| 2 Does not include parameter sets producing C:R > 10. | ||||||
CRsummProd <- CRcomparison_tables %>%
ungroup() %>%
filter(., Function == "Trophic Productivity") %>%
dplyr::mutate(., Compartment = fct_recode(Compartment, Consumers = "PROD_C1rr",
Consumers = "PROD_C2rr", Consumers = "PROD_Crr", Producers = "PROD_P1rr",
Producers = "PROD_P2rr", Producers = "PROD_Prr")) %>%
dplyr::mutate(., Compartment = fct_relevel(Compartment, "Producers",
"Consumers")) %>%
pivot_wider(., names_from = c(CR, Fertility), values_from = median) %>%
gt(., groupname_col = "Scale") %>%
fmt_number(., columns = 4:9, decimals = 2) %>%
tab_spanner(label = md("**Along-gradient**"), columns = c("all_Along",
"above_1_Along", "below_1_Along")) %>%
tab_spanner(label = md("**Against-gradient**"), columns = c("all_Against",
"above_1_Against", "below_1_Against")) %>%
cols_label(Compartment = md("**Compartment**"), all_Against = md("*All C:R*"),
above_1_Against = md("*C:R > 1*"), below_1_Against = md("*C:R < 1*"),
all_Along = md("*All C:R*"), above_1_Along = md("*C:R > 1*"), below_1_Along = md("*C:R < 1*")) %>%
cols_hide("Function") %>%
tab_footnote(footnote = "Includes parameter sets producing C:R > 10.",
locations = cells_column_labels(columns = c("all_Against", "all_Along"))) %>%
tab_footnote(footnote = "Does not include parameter sets producing C:R > 10.",
locations = cells_column_labels(columns = c("above_1_Against",
"above_1_Along"))) %>%
tab_header(title = md("**Median LRR values of trophic compartment productivity**"),
subtitle = md("LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing _C:R_ > 1, or only parameter sets producing _C:R_ < 1.")) %>%
tab_style(style = cell_text(align = "left", size = 18), locations = cells_title())
CRsummProd
| Median LRR values of trophic compartment productivity | ||||||
| LRR values are grouped by fertility scenario and by whether they were calculated using all parameter sets, only parameter sets producing C:R > 1, or only parameter sets producing C:R < 1. | ||||||
| Compartment | Against-gradient | Along-gradient | ||||
|---|---|---|---|---|---|---|
| All C:R1 | C:R > 12 | C:R < 1 | All C:R1 | C:R > 12 | C:R < 1 | |
| Ecosystem 1 | ||||||
| Consumers | −0.72 | −0.71 | −0.72 | 0.26 | 0.26 | 0.26 |
| Producers | −0.16 | −0.55 | −0.15 | 0.12 | 0.20 | 0.08 |
| Ecosystem 2 | ||||||
| Consumers | 0.22 | 0.20 | 0.22 | −0.47 | −0.28 | −0.57 |
| Producers | 0.18 | 0.17 | 0.16 | −0.34 | −0.27 | −0.34 |
| Meta-ecosystem | ||||||
| Consumers | 0.03 | −0.07 | −0.14 | −0.04 | 0.06 | −0.05 |
| Producers | 0.03 | 0.01 | 0.03 | −0.04 | −0.02 | −0.05 |
| 1 Includes parameter sets producing C:R > 10. | ||||||
| 2 Does not include parameter sets producing C:R > 10. | ||||||
Here, we perform a Global Sensitivity Analysis (GSA; sensu Harper, Stella, and Fremier (2011)) on the parameter sets that used to produce our results—i.e., those resulting in stable equilibria with biological meaning. Compared to local sensitivity analyses, whereby each parameter is varied over a range of values one at the time while all others are kept constant, GSA allows for varying all parameters at once. In turn, this has several advantages, namely (i) it removes the need to elect a constant value at which all parameters are kept while varying the one under analysis and (ii) it accounts for parameter interactions and is thus particularly suited for complex mathematical models with multiple parameters and non-linear relationships.
We use function randomForest() from package randomForest to run our sensitivity
analysis. Rather than running a single analysis on the whole dataset, we run two
separate ones—one each for the donor ecosystem 1 and the recipient ecosystem 2.
This allows us to analyse parameter importance in each consumer movement
scenario—gradient-neutral, along-, and against-gradient—and across
the four different biomass distribution conditions identified above: (i)
bottom-heavy pyramid in the donor ecosystem 1, (ii) bottom-heavy pyramid in the
recipient ecosystem 2, (iii) top-heavy pyramid in the donor ecosystem 1, (iv)
top-heavy pyramid in the recipient ecosystem 2
First, we load the ecoCR object, and create two additional variables:
E1CRabove1 and E2CRabove1. Both these are categorical variables, with Yes
or No values, and capture whether the local biomass distribution in the donor
ecosystem 1 or in the recipient ecosystem 2 conforms to a classic bottom-heavy
pyramid or to an inverted, top-heavy pyramid.
Afterwards, we store the names of the six state variables in our model into a
new object stateV that we will use later in the code to isolate each state
variable in turn.
# load the data, then create additional categorical variables for
# plotting
ecoCR <- readRDS(file = "../Results/ecoCR.rds") %>%
mutate(., E1CRabove1 = as_factor(ifelse(Eco1 > 1, "Y", "N")), E2CRabove1 = as_factor(ifelse(Eco2 >
1, "Y", "N")), Fertility = as_factor(Fertility))
# define a vector to use in the for loop to go through each state
# variable in the model
stateV <- c("N1", "P1", "C1", "N2", "P2", "C2")
We will run two Global Sensitivity Analyses on the whole ecoCR object, based
on whether we take into account the biomass distribution in the donor ecosystem
1 or in the recipient ecosystem 2. This is necessary as the same parameter set
can give rise to an inverted, top-heavy biomass pyramid in either local
ecosystem but to a classic, bottom-heavy one in the other. By focusing on one
ecosystem at the time in terms of biomass distribution, we can account for all
four scenarios highlighted above.
Starting with the donor ecosystem 1, we use a for loop to run out GSA on our
ecoCR object, one state variable at the time. For each run of the loop, we
isolate the relevant state variable, generate the random forest trees, calculate
importance and relative importance, and store all this in a new object,
sensiResE1.
Results for this analysis are shown in Figures 27 and 29.
# allocate an empty dataframe to contain the results of each
# sensitivity analysis
sensiResE1 <- NULL
# This loop runs the sensitivity analysis
for (i in 1:length(stateV)) {
# wrangle the dataframe to isolate the columns we want to work
# with Note that we are only include E1CRabove1 because here we
# focus on biomass distribution (i.e., C:R value) in the donor
# ecosystem 1
sensiA <- select(ecoCR, stateV[i], I1:e2, Fertility, E1CRabove1)
# rename the state variable column to a neutral name
colnames(sensiA)[1] <- "Compartment"
# here I group by Fertility (= movement scenario) and whether the
# CR in Ecosystem 1 (donor) is > 1, so I can later use these two
# pieces of information to produce multiple graphs from the same
# df (which reduces the number of runs for the for loop) The last
# function call in the pine runs the randomForest model to
# estimate
sA_temp <- sensiA %>%
group_by(Fertility, E1CRabove1) %>%
select(., Compartment:e2) %>%
do(Cmove = randomForest(Compartment ~ . - Fertility - E1CRabove1,
data = ., proximity = TRUE, importance = TRUE)) # run this and see
# allocate empty dataframe to store results from next for loop
imp_temp <- NULL
# this loop gets the results for each movement scenario
for (j in seq(1, 6, 1)) {
# this pipe extracts the importance values of each parameter
# and stores them in the df created just above
temp <- as_tibble(importance(sA_temp[[3]][[j]]), rownames = NA) %>%
rownames_to_column(., var = "Parameter") %>%
mutate(., Movement = as_factor(sA_temp[[1]][[j]]), E1CRabove1 = as_factor(sA_temp[[2]][[j]]))
imp_temp <- rbind(imp_temp, temp)
# browser() print(j)
}
# now calculate the ranked importance value as the ratio of each
# parameter's importance values to the sum of all parameters'
# importance values, store it in a temporary object to which we
# add grouping and categorical variables that are useful in
# plotting the results
ranked_imp_temp <- tibble(imp_temp[, 3]/sum(imp_temp[, 3])) %>%
mutate(., Parameter = as_factor(imp_temp$Parameter), Movement = as_factor(imp_temp$Movement),
E1CRabove1 = as_factor(imp_temp$E1CRabove1), SV = as_factor(stateV[i]))
# store the results from each loop run in a separate df so the
# loop can run again without loosing results
sensiResE1 <- bind_rows(sensiResE1, ranked_imp_temp)
# browser()
print(i)
}
[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
The GSA taking into account the biomass distribution of ecosystem 2 runs along
the same lines as the one for the donor ecosystem 1, with the results being
stored in object sensiResE2.
Results for this analysis are shown in Figures 28 and 30.
# allocate an empty dataframe to contain the results of each
# sensitivity analysis
sensiResE2 <- NULL
# This loop runs the sensitivity analysis
for (i in 1:length(stateV)) {
# wrangle the dataframe to isolate the columns we want to work
# with Note that we are only include E2CRabove1 because here we
# focus on biomass distribution (i.e., C:R value) in the
# recipient ecosystem 2
sensiA <- select(ecoCR, stateV[i], I1:e2, Fertility, E2CRabove1)
# rename the state variable column to a neutral name
colnames(sensiA)[1] <- "Compartment"
# here I group by Fertility (= movement scenario) and whether the
# CR in Ecosystem 2 (recipient) is > 1, so I can later use these
# two pieces of information to produce multiple graphs from the
# same df (which reduces the number of runs for the for loop) The
# last function call in the pine runs the randomForest model to
# estimate
sA_temp <- sensiA %>%
group_by(Fertility, E2CRabove1) %>%
select(., Compartment:e2) %>%
do(Cmove = randomForest(Compartment ~ . - Fertility - E2CRabove1,
data = ., proximity = TRUE, importance = TRUE)) # run this and see
# allocate empty dataframe to store results from next for loop
imp_temp <- NULL
# this loop gets the results for each movement scenario
for (j in seq(1, 6, 1)) {
# this pipe extracts the importance values of each parameter
# and stores them in the df created just above
temp <- as_tibble(importance(sA_temp[[3]][[j]]), rownames = NA) %>%
rownames_to_column(., var = "Parameter") %>%
mutate(., Movement = as_factor(sA_temp[[1]][[j]]), E2CRabove1 = as_factor(sA_temp[[2]][[j]]))
imp_temp <- rbind(imp_temp, temp)
# browser() print(j)
}
# now calculate the ranked importance value as the ratio of each
# parameter's importance values to the sum of all parameters'
# importance values, store it in a temporary object to which we
# add grouping and categorical variables that are useful in
# plotting the results
ranked_imp_temp <- tibble(imp_temp[, 3]/sum(imp_temp[, 3])) %>%
mutate(., Parameter = as_factor(imp_temp$Parameter), Movement = as_factor(imp_temp$Movement),
E2CRabove1 = as_factor(imp_temp$E2CRabove1), SV = as_factor(stateV[i]))
# store the results from each loop run in a separate df so the
# loop can run again without loosing results
sensiResE2 <- bind_rows(sensiResE2, ranked_imp_temp)
# browser()
print(i)
}
[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
As shown in the figures below (Figures @ref(fig: gsa-results-CR-below-1-eco-1) to 30), the shape of the biomass pyramid in either the donor or recipient ecosystem does not appear to influence the relative importance of the parameters in the model. For this reason, we run a final GSA on the whole model. Results from this model run are shown in Figure 26.
# allocate an empty dataframe to contain the results of each
# sensitivity analysis
sensiRes <- NULL
# This loop runs the sensitivity analysis
for (i in 1:length(stateV)) {
# wrangle the dataframe to isolate the columns we want to work
# with Note that we are only include E2CRabove1 because here we
# focus on biomass distribution (i.e., C:R value) in the
# recipient ecosystem 2
sensiA <- select(ecoCR, stateV[i], I1:e2, Fertility)
# rename the state variable column to a neutral name
colnames(sensiA)[1] <- "Compartment"
# here I group by Fertility (= movement scenario) and whether the
# CR in Ecosystem 2 (recipient) is > 1, so I can later use these
# two pieces of information to produce multiple graphs from the
# same df (which reduces the number of runs for the for loop) The
# last function call in the pine runs the randomForest model to
# estimate
sA_temp <- sensiA %>%
group_by(Fertility) %>%
select(., Compartment:e2) %>%
do(Cmove = randomForest(Compartment ~ . - Fertility, data = .,
proximity = TRUE, importance = TRUE)) # run this and see
# allocate empty dataframe to store results from next for loop
imp_temp <- NULL
# this loop gets the results for each movement scenario
for (j in seq(1, 3, 1)) {
# this pipe extracts the importance values of each parameter
# and stores them in the df created just above
temp <- as_tibble(importance(sA_temp[[2]][[j]]), rownames = NA) %>%
rownames_to_column(., var = "Parameter") %>%
mutate(., Movement = as_factor(sA_temp[[1]][[j]]))
imp_temp <- rbind(imp_temp, temp)
# browser() print(j)
}
# now calculate the ranked importance value as the ratio of each
# parameter's importance values to the sum of all parameters'
# importance values, store it in a temporary object to which we
# add grouping and categorical variables that are useful in
# plotting the results
ranked_imp_temp <- tibble(imp_temp[, 3]/sum(imp_temp[, 3])) %>%
mutate(., Parameter = as_factor(imp_temp$Parameter), Movement = as_factor(imp_temp$Movement),
SV = as_factor(stateV[i]))
# store the results from each loop run in a separate df so the
# loop can run again without loosing results
sensiRes <- bind_rows(sensiRes, ranked_imp_temp)
# browser()
print(i)
}
[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6
Now that we have a relative importance value for each parameter in each movement scenario, in all four conditions of biomass distribution in the two local ecosystems (see above, Changes to biomass distribution in local and meta-ecosystem), we produce graphs to visualize these results. These graphs will be useful in identifying the more important parameters in determining the stock value of each trophic compartment in the model for each combination of movement scenario and biomass distribution conditions.
Figure 26 presents the results from the GSA run on the whole system, which do not differ from those of two GSA runs on the two local ecosystems (see above; Whole System). Figures 27 to 30 below present the results from the local-ecosystem GSAs, for transparancy and completeness.
# plotting pipe first, rename the variable to use as x-axis, then
# filter by C:R > 1 in the recipient ecosystem 2, then create a new
# variable to store information about which ecosystem do parameters
# belong to. Then, turn it into a factor, reorder the levels using
# fct_relevel(), and reorder the levels of factor Parameter and their
# associated relative importance values according to the Scale
# factor. Finally, group by Scale and plot
GSAwhole <- sensiRes %>%
rename(., RelImpI = IncNodePurity) %>%
mutate(., Scale = ifelse(Parameter == "u1" | Parameter == "h1" | Parameter ==
"a1" | Parameter == "e1" | Parameter == "d1" | Parameter == "I1",
"Donor", ifelse(Parameter == "g" | Parameter == "m" | Parameter ==
"c", "Dispersers' Pool", ifelse(Parameter == "l", "Other",
"Recipient"))), Scale = as_factor(Scale), Scale = fct_relevel(Scale,
"Donor", "Recipient", "Dispersers' Pool", "Other"), Parameter = fct_reorder2(.f = Parameter,
.x = RelImpI, .y = Scale)) %>%
group_by(., Scale) %>%
ggplot(., aes(x = RelImpI, y = Parameter)) + geom_col(aes(fill = Scale),
col = "black", linewidth = 0.25) + scale_fill_met_d(name = "Egypt",
direction = 1) + xlab("Relative Importance Index") + ylab(" ") + facet_grid2(Movement ~
SV) + theme(legend.position = "top", legend.title = element_blank(),
text = element_text(size = 20), panel.spacing = unit(0.5, "lines"))
GSAwhole + scale_y_discrete("Parameter", labels = c("l", "g", "m", "c",
expression(I[2]), expression(u[2]), expression(a[2]), expression(h[2]),
expression(d[2]), expression(e[2]), expression(I[2]), expression(u[1]),
expression(a[1]), expression(h[1]), expression(d[1]), expression(e[1])))
Figure 26: At the whole system level, considering together all cases of biomass distribution in the two ecosystems, the global sensitivity analysis again shows evidence for the direct and indirect effects of functional response-related parameters across all movement scenarios and in both ecosystems. In particular, when consumers move against-gradient (bottom row), notice the importance of parameters e2 and d2 to shape not only the dynamics of trophic compartment in the recipient ecosystem 2 but also those of the primary producers in ecosystem 1. Note that this figure is equivalent to Figure A.7 in the Appendices to the main text, but with x-axes tailored to each panel’s range of relative importance values.
# save statement
ggsave("../Results/SA_RelImpI_WS.pdf", device = "pdf", dpi = 300, width = 12,
height = 10)
# plotting pipe first, rename the variable to use as x-axis, then
# filter by C:R < 1 in the donor ecosystem 1, then create a new
# variable to store information about which ecosystem do parameters
# belong to. Then, turn it into a factor, reorder the levels using
# fct_relevel(), and reorder the levels of factor Parameter and their
# associated relative importance values according to the Scale
# factor. Finally, group by Scale and plot
sensiResE1 %>%
rename(., RelImpI = IncNodePurity) %>%
filter(., E1CRabove1 == "N") %>%
mutate(., Scale = ifelse(Parameter == "u1" | Parameter == "h1" | Parameter ==
"a1" | Parameter == "e1" | Parameter == "d1" | Parameter == "I1",
"Donor", ifelse(Parameter == "g" | Parameter == "m" | Parameter ==
"c", "Dispersers' Pool", ifelse(Parameter == "l", "Other",
"Recipient"))), Scale = as_factor(Scale), Scale = fct_relevel(Scale,
"Donor", "Recipient", "Dispersers' Pool", "Other"), Parameter = fct_reorder2(.f = Parameter,
.x = RelImpI, .y = Scale)) %>%
group_by(., Scale) %>%
ggplot(., aes(x = RelImpI, y = Parameter)) + geom_col(aes(fill = Scale),
col = "black", linewidth = 0.25) + scale_fill_met_d(name = "Egypt",
direction = 1) + xlab("Relative Importance Index") + ylab(" ") + facet_grid2(Movement ~
SV) + theme(legend.position = "top", legend.title = element_blank(),
text = element_text(size = 20), panel.spacing = unit(1.5, "lines")) +
scale_y_discrete("Parameter", labels = c("l", "g", "m", "c", expression(I[2]),
expression(u[2]), expression(a[2]), expression(h[2]), expression(d[2]),
expression(e[2]), expression(I[2]), expression(u[1]), expression(a[1]),
expression(h[1]), expression(d[1]), expression(e[1])))
Figure 27: The global sensitivity analysis (GSA) performed on the subset of parameter sets that produces bottom-heavy biomass pyramids in the donor ecosystem provides key insights in the direct and indirect pathways through which higher trophic compartment influence lower ones. The high relative importance of parameter u1, involved in the functional response of Primary Producers, is a clear example of a direct effect shaping the accumulation of Nutrients in the donor ecosystem 1 in all three movement scenario (first column). Likewise, the high relative importance of parameter e2, the Consumer efficiency in the recipient ecosystem 2, in shaping Nutrient stocks in the recipient ecosystem 2 (fourth column, bottom row) highlights the strength of the indirect pathway that allows the highest trophic level in our system to influence the lower one. Note that e2 is also the most importance parameter in determining biomass accumulation of Consumers in the recipient ecosystem 2 in this scenario, meaning that it has both direct and indirect effects on the system’s dynamics.
# save statement
ggsave("../Results/SA_RelImpI_Eco1CRbelow1.pdf", device = "pdf", dpi = 300,
width = 12, height = 10)
# plotting pipe first, rename the variable to use as x-axis, then
# filter by C:R < 1 in the recipient ecosystem 2, then create a new
# variable to store information about which ecosystem do parameters
# belong to. Then, turn it into a factor, reorder the levels using
# fct_relevel(), and reorder the levels of factor Parameter and their
# associated relative importance values according to the Scale
# factor. Finally, group by Scale and plot
sensiResE2 %>%
rename(., RelImpI = IncNodePurity) %>%
filter(., E2CRabove1 == "N") %>%
mutate(., Scale = ifelse(Parameter == "u1" | Parameter == "h1" | Parameter ==
"a1" | Parameter == "e1" | Parameter == "d1" | Parameter == "I1",
"Donor", ifelse(Parameter == "g" | Parameter == "m" | Parameter ==
"c", "Dispersers' Pool", ifelse(Parameter == "l", "Other",
"Recipient"))), Scale = as_factor(Scale), Scale = fct_relevel(Scale,
"Donor", "Recipient", "Dispersers' Pool", "Other"), Parameter = fct_reorder2(.f = Parameter,
.x = RelImpI, .y = Scale)) %>%
group_by(., Scale) %>%
ggplot(., aes(x = RelImpI, y = Parameter)) + geom_col(aes(fill = Scale),
col = "black", linewidth = 0.25) + scale_fill_met_d(name = "Egypt",
direction = 1) + xlab("Relative Importance Index") + ylab(" ") + facet_grid2(Movement ~
SV) + theme(legend.position = "top", legend.title = element_blank(),
text = element_text(size = 20), panel.spacing = unit(1.5, "lines")) +
scale_y_discrete("Parameter", labels = c("l", "g", "m", "c", expression(I[2]),
expression(u[2]), expression(a[2]), expression(h[2]), expression(d[2]),
expression(e[2]), expression(I[2]), expression(u[1]), expression(a[1]),
expression(h[1]), expression(d[1]), expression(e[1])))
Figure 28: When the recipient ecosystem 2 has a bottom-heavy biomass pyramid, we see evidence of direct and indirect effects of high trophic compartments on lower ones arising from our global sensitivity analysis. Of particular note are the direct and indirect effects of parameter e1 on the biomass accumulation of Consumers and Nutrient in the donor ecosystem 1, respectively, when consumer move along-gradient.
# save statement
ggsave("../Results/SA_RelImpI_Eco2CRbelow1.pdf", device = "pdf", dpi = 300,
width = 12, height = 10)
# plotting pipe first, rename the variable to use as x-axis, then
# filter by C:R > 1 in the donor ecosystem 1, then create a new
# variable to store information about which ecosystem do parameters
# belong to. Then, turn it into a factor, reorder the levels using
# fct_relevel(), and reorder the levels of factor Parameter and their
# associated relative importance values according to the Scale
# factor. Finally, group by Scale and plot
sensiResE1 %>%
rename(., RelImpI = IncNodePurity) %>%
filter(., E1CRabove1 == "Y") %>%
mutate(., Scale = ifelse(Parameter == "u1" | Parameter == "h1" | Parameter ==
"a1" | Parameter == "e1" | Parameter == "d1" | Parameter == "I1",
"Donor", ifelse(Parameter == "g" | Parameter == "m" | Parameter ==
"c", "Dispersers' Pool", ifelse(Parameter == "l", "Other",
"Recipient"))), Scale = as_factor(Scale), Scale = fct_relevel(Scale,
"Donor", "Recipient", "Dispersers' Pool", "Other"), Parameter = fct_reorder2(.f = Parameter,
.x = RelImpI, .y = Scale)) %>%
group_by(., Scale) %>%
ggplot(., aes(x = RelImpI, y = Parameter)) + geom_col(aes(fill = Scale),
col = "black", linewidth = 0.25) + scale_fill_met_d(name = "Egypt",
direction = 1) + xlab("Relative Importance Index") + ylab(" ") + facet_grid2(Movement ~
SV) + theme(legend.position = "top", legend.title = element_blank(),
text = element_text(size = 20), panel.spacing = unit(1.5, "lines")) +
scale_y_discrete("Parameter", labels = c("l", "g", "m", "c", expression(I[2]),
expression(u[2]), expression(a[2]), expression(h[2]), expression(d[2]),
expression(e[2]), expression(I[2]), expression(u[1]), expression(a[1]),
expression(h[1]), expression(d[1]), expression(e[1])))
Figure 29: When the donor ecosystem 1 has a top-heavy biomass pyramid, we note again strong direct effects of functional response-related parameters on the dynamics of Consumers and Nutrient in this ecosystem, when consumers move either in a gradient-neutral or along-gradient fashion. Note that weak indirect effects of parameter e1 on Nutrient stock accumulation in the donor ecosystem 1 are still visible when consumers move along-gradient (first column, middle row). In the against-gradient movement scenario, we note again some evidence of indirect effects of parameter e2 on the dynamics of the Nutrient compartment in the recipient ecosystem 2 (fourth column, bottom row).
# save statement
ggsave("../Results/SA_RelImpI_Eco1CRabove1.pdf", device = "pdf", dpi = 300,
width = 12, height = 10)
# plotting pipe first, rename the variable to use as x-axis, then
# filter by C:R > 1 in the recipient ecosystem 2, then create a new
# variable to store information about which ecosystem do parameters
# belong to. Then, turn it into a factor, reorder the levels using
# fct_relevel(), and reorder the levels of factor Parameter and their
# associated relative importance values according to the Scale
# factor. Finally, group by Scale and plot
sensiResE2 %>%
rename(., RelImpI = IncNodePurity) %>%
filter(., E2CRabove1 == "Y") %>%
mutate(., Scale = ifelse(Parameter == "u1" | Parameter == "h1" | Parameter ==
"a1" | Parameter == "e1" | Parameter == "d1" | Parameter == "I1",
"Donor", ifelse(Parameter == "g" | Parameter == "m" | Parameter ==
"c", "Dispersers' Pool", ifelse(Parameter == "l", "Other",
"Recipient"))), Scale = as_factor(Scale), Scale = fct_relevel(Scale,
"Donor", "Recipient", "Dispersers' Pool", "Other"), Parameter = fct_reorder2(.f = Parameter,
.x = RelImpI, .y = Scale)) %>%
group_by(., Scale) %>%
ggplot(., aes(x = RelImpI, y = Parameter)) + geom_col(aes(fill = Scale),
col = "black", linewidth = 0.25) + scale_fill_met_d(name = "Egypt",
direction = 1) + xlab("Relative Importance Index") + ylab(" ") + facet_grid2(Movement ~
SV) + theme(legend.position = "top", legend.title = element_blank(),
text = element_text(size = 20), panel.spacing = unit(1.5, "lines")) +
scale_y_discrete("Parameter", labels = c("l", "g", "m", "c", expression(I[2]),
expression(u[2]), expression(a[2]), expression(h[2]), expression(d[2]),
expression(e[2]), expression(I[2]), expression(u[1]), expression(a[1]),
expression(h[1]), expression(d[1]), expression(e[1])))
Figure 30: When the recipient ecosystem 2 has a top-heavy biomass pyramid, the high relative importance of functional response-related parameters in both the dynamics of the Consumer and Nutrient compartments of the recipient ecosystem 2 points to both direct and indirect effects being exerted by Consumers in the against-gradient and gradient-neutral movement scenarios (forth and sixth columns, top and bottom rows).
# save statement
ggsave("../Results/SA_RelImpI_Eco2CRabove1.pdf", device = "pdf", dpi = 300,
width = 12, height = 10)
This notebook was compiled in 37 minutes.
We are grateful to Dr. Anne McLeod for her help in developing the stability and sensitivity analyses code, and for suggesting useful reference on both topics. We thank Dr. Robert Buchkowski for his insights in designing and evaluating the stability analyses. Finally, we thank Dr. A. Z. Andis Arietta for his help in developing the code for the biomass pyramid plots.
MR and SJL conceptualized the model; MR, SJL, OJS, EVW, YFW, and TRH designed the study; MR wrote the code and run the analyses; MR and SJL run and proofread the code; MR, SJL, OJS interpreted the results; all authors contributed critically to the manuscript.
If you see mistakes or want to suggest changes, please create an issue on the source repository.
Text and figures are licensed under Creative Commons Attribution CC BY-SA 4.0. Source code is available at https://github.com/matteorizzuto/AnimalVectorNutrientFlow, unless otherwise noted. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".